[线性代数] 矩阵求逆 不知道哪一步错了

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

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

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