矩阵求逆 不知道哪一步错了

kerb

下面的程序是一位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: [url=mailto:CHUGUOZHEN@ALIYUN.COM]CHUGUOZHEN@ALIYUN.COM[/url]

     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