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

fcode

[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