Fortran Coder
标题: 矩阵求逆 不知道哪一步错了 [打印本页]
作者: gdx1958861161 时间: 2015-3-17 10:28
标题: 矩阵求逆 不知道哪一步错了
[Fortran] 纯文本查看 复制代码
PROGRAM QIUNI
PARAMETER(NN=3)
REAL*8 ff(NN,NN),iff(NN,NN)
ff(1,1)=2; ff(1,2)=2; ff(1,3)=1
ff(2,1)=3; ff(2,2)=1; ff(2,3)=5
ff(3,1)=3; ff(3,2)=2; ff(3,3)=3
CALL QN(ff,iff,NN)
do i=1,3
do j=1,3
print*,iff(i,j)
enddo
enddo
END
!***********************************************************************
!* 子程序1：调用高斯法求逆矩阵 *
!***********************************************************************
!**********************************************************************
SUBROUTINE QN(A1,AN,N)
REAL*8 A1(N,N),A2(N,N),AN(N,N),B(N),X(N)
AN=0D0
B=0D0
DO 10 J=1,N
A2=A1
B(J)=1
CALL GS(A2,B,N,X)
AN(1:N,J)=X(1:N)
10 CONTINUE
END
!***********************************************************************
!* 子程序2：列主元素Gauss消去法 *
!***********************************************************************
SUBROUTINE GS(A,B,N,X)
REAL*8 A(N,N),AMAX(1,N),B(N),X(N),M,BB,MAX,SUM
DO 40 K=1,N-1
MAX=DABS(A(K,K))
L=K
DO 10 I=K,N
IF(DABS(A(I,K))>=MAX) THEN
MAX=DABS(A(I,K))
L=I
ENDIF
10 CONTINUE
AMAX(1,K:N)=A(L,K:N)
A(L,K:N)=A(K,K:N)
A(K,K:N)=AMAX(1,K:N)
BB=B(L)
B(L)=B(K)
B(K)=BB
DO 30 I=K+1,N
M=A(I,K)/A(K,K)
DO 20 J=K,N
A(I,J)=A(I,J)-M*A(K,J)
20 CONTINUE
B(I)=B(I)-M*B(K)
30 CONTINUE
40 CONTINUE
X(N)=B(N)/A(N,N)
DO 60 K=N-1,1,-1
SUM=0
DO 50 J=K+1,N
SUM=SUM+A(K,J)*X(J)
50 CONTINUE
X(K)=(B(K)-SUM)/A(K,K)
60 CONTINUE
END
作者: gdx1958861161 时间: 2015-3-17 10:29
能不能提供一个可以正确计算的矩阵求逆程序
作者: fcode 时间: 2015-3-17 12:39
[Fortran] 纯文本查看 复制代码
Program www_fcode_cn
Implicit None
Integer , Parameter :: nn=3
Real(Kind=8) :: ff(nn, nn), iff(nn, nn)
integer :: j
ff(:,:) = reshape( [2,3,3,2,1,2,1,5,3] , [nn,nn] )
Call qn(ff, iff, nn)
Do j = 1 , nn
print*,iff(:,j)
End Do
Write (*, *) matmul(ff, iff) !// 验证
End Program www_fcode_cn
!* 子程序1：调用高斯法求逆矩阵
Subroutine qn(a1, an, n)
Implicit None
Integer :: n
Real(Kind=8) :: a1(n, n), a2(n, n), an(n, n), b(n), x(n)
integer :: j
an = 0.0d0
Do j = 1, n
a2 = a1
b = 0.D0 !// 此处放入循环体内
b(j) = 1.0d0
Call gs(a2, b, n, x)
an(1:n, j) = x(1:n)
End Do
End Subroutine qn
!* 子程序2：列主元素Gauss消去法
Subroutine gs(a, b, n, x)
Implicit None
Integer :: N
Real(Kind=8) ::a(n, n), amax(1, n), b(n), x(n), m, bb, max, sum
integer :: i , j , k , l
Do k = 1, n - 1
max = dabs(a(k,k))
l = k
Do i = k, n
If (dabs(a(i,k))>=max) Then
max = dabs(a(i,k))
l = i
End If
End Do
amax(1, k:n) = a(l, k:n)
a(l, k:n) = a(k, k:n)
a(k, k:n) = amax(1, k:n)
bb = b(l)
b(l) = b(k)
b(k) = bb
Do i = k + 1, n
m = a(i, k)/a(k, k)
a(i, k:n) = a(i, k:n) - m*a(k, k:n)
b(i) = b(i) - m*b(k)
End Do
End Do
x(n) = b(n)/a(n, n)
Do k = n - 1, 1, -1
sum = 0.0d0
Do j = k + 1, n
sum = sum + a(k, j)*x(j)
End Do
x(k) = (b(k)-sum)/a(k, k)
End Do
End Subroutine gs
作者: kerb 时间: 2015-3-21 22:43
本帖最后由 kerb 于 2016-8-12 18:41 编辑
下面的程序是一位50多岁的老人写的，刚学fortran不久后的作品，可以参考一下
[Fortran] 纯文本查看 复制代码
PROGRAM MAIN
IMPLICIT NONE
INTEGER ( KIND = 8 ) :: N = 1000
REAL(8) ELAPSED_TIME
REAL(8), DIMENSION( 1000, 1000 ) :: A
REAL(8), DIMENSION( 1000 ) :: B
REAL(8), DIMENSION( 1000 ) :: X
INTEGER ( KIND = 8 ) I
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(B)
CALL RANDOM_NUMBER(A)
WRITE ( UNIT = 6, FMT = * ) "EXAMPLE OF GAUSSIAN ELIMINATION AND BACK SUBSTITUTION: A * X = B"
WRITE ( UNIT = 6, FMT = * ) "SAVE MATRIX TO FILE"
OPEN( UNIT = 1, FILE='MAT.TXT')
WRITE ( UNIT = 1, FMT = * ) A
CLOSE ( 1 )
WRITE ( UNIT = 6, FMT = * ) "SAVE RHS TO FILE"
OPEN ( UNIT = 2, FILE = 'RHS.TXT' )
WRITE ( UNIT = 2, FMT = * ) B
CLOSE ( 2 )
WRITE ( UNIT = 6, FMT = * ) "SOLVING......"
!ELAPSED_TIME = TIMEF( )
!CALL GAUSSJ ( A, N, B, X )
OPEN( UNIT = 3, FILE='SOL.TXT')
!ELAPSED_TIME = TIMEF( )
!WRITE ( UNIT = 6, FMT = * ) "ELAPSED TIME: ", ELAPSED_TIME, "SECONDS"
WRITE ( UNIT = 6, FMT = * ) "SAVE SOLUTION TO FILE"
WRITE ( UNIT = 3, FMT = * ) X
CLOSE ( 3 )
WRITE ( UNIT = 6, FMT = * ) "SOLVING HAS COMPLETED, PRESS ENTER TO RETURN"
PAUSE
END PROGRAM
SUBROUTINE GAUSSJ ( A, N, B, X )
!AUTHOR: G. Z. CHU
!EMAIL_TO: CHUGUOZHEN@ALIYUN.COM
IMPLICIT NONE
INTEGER ( KIND = 8 ) :: N
REAL(8) :: A ( N, N )
REAL(8) :: B ( N )
REAL(8) :: X ( N )
INTEGER ( KIND = 8 ) :: IPIV ( N )
INTEGER ( KIND = 8 ) :: IT, JT, KT, I, J, K, M, INDEX
REAL(8) :: TEMP2, TEMP3
REAL(8) :: MAX
DO I = 1, N
IPIV ( I ) = I
END DO
DO I = 1, N
IT = IPIV ( I )
MAX = DABS ( A ( I, IT ) )
INDEX = I
DO J = I + 1, N
JT = IPIV ( I )
IF ( MAX < DABS ( A ( I, JT ) ) ) THEN
MAX = DABS ( A ( I, JT ) )
INDEX = J
END IF
END DO
IF ( INDEX /= IT ) THEN
KT = IPIV ( I )
IPIV ( I ) = IPIV ( INDEX )
IPIV ( INDEX ) = KT
ENDIF
IT = IPIV ( I )
TEMP2 = 1.D0 / A ( I, IT )
DO K = I + 1, N
KT = IPIV ( K )
!IF ( A ( I, KT ) /= 0.D0 ) THEN !IF "A" IS A SPARSE MATRIX
TEMP3 = A ( I, KT ) * TEMP2
DO M = I, N
!IF ( A ( M, IT ) /= 0.D0 ) THEN !IF "A" IS A SPARSE MATRIX
A ( M, KT ) = A ( M, KT ) - A ( M, IT ) * TEMP3
!ENDIF
END DO
B ( KT ) = B ( KT ) - B ( IT ) * TEMP3
!ENDIF
END DO
END DO
!BACK-SUBSTITUTION
DO I = N, 1, -1
TEMP3 = 0
IT = IPIV ( I )
TEMP2 = 1.D0 / A ( I, IT )
DO K = I + 1, N
!IF ( A ( K, IT ) /= 0.D0 ) THEN !IF "A" IS A SPARSE MATRIX
TEMP3 = TEMP3 + A ( K, IT ) * X ( K )
!END IF
END DO
X ( I ) = ( B ( IT ) -TEMP3 ) * TEMP2
END DO
END SUBROUTINE
作者: cherish_J 时间: 2019-3-29 09:44
你好！我调用这个子程序后，出现堆栈溢出的错误。计算的矩阵维数是1024*1024的，请问应该怎么解决呢？我搜了堆栈溢出错误的解决方法，还是不知道错在哪儿
作者: 千手斑 时间: 2019-4-6 23:41
您好，请问我现在想在列主元消去法中用openmp实现并行，但是看其中的循环嵌套有点复杂，这好实现吗，我对openmp理解的还不好。
作者: aoeo2jam 时间: 2019-5-26 15:48
最好用全主元消元，高斯消元只适用于小规模的线性代数方程组，大规模情况下舍入误差累计易导致计算失败
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) |
Powered by Discuz! X3.2 |