三种方法用Fortran求逆矩阵

三种方法用Fortran求四阶矩阵的逆矩阵

数值计算

C     refer to https://fortranwiki.org/fortran/show/Matrix+inversion
      SUBROUTINE MATINV(A,B)
      DIMENSION  A(4,4),B(4,4)
      d = 1/
     &  (A(1,1)*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))
     &  +A(2,3)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(2,4)*
     &  (A(3,2)*A(4,3)-A(3,3)*A(4,2)))
     & - A(1,2)*(A(2,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))
     &  +A(2,3)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*
     &  (A(3,1)*A(4,3)-A(3,3)*A(4,1)))
     & + A(1,3)*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))
     &  +A(2,2)*(A(3,4)*A(4,1)-A(3,1)*A(4,4))+A(2,4)*
     &  (A(3,1)*A(4,2)-A(3,2)*A(4,1)))
     & - A(1,4)*(A(2,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))
     &  +A(2,2)*(A(3,3)*A(4,1)-A(3,1)*A(4,3))+A(2,3)*
     &  (A(3,1)*A(4,2)-A(3,2)*A(4,1))))
      B(1,1) = d*(A(2,2)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(2,3)*(A(3,4)
     &   *A(4,2)-A(3,2)*A(4,4))+A(2,4)*(A(3,2)*A(4,3)-A(3,3)*A(4,2)))
      B(2,1) = d*(A(2,1)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(2,3)*(A(3,1)
     &   *A(4,4)-A(3,4)*A(4,1))+A(2,4)*(A(3,3)*A(4,1)-A(3,1)*A(4,3)))
      B(3,1) = d*(A(2,1)*(A(3,2)*A(4,4)-A(3,4)*A(4,2))+A(2,2)*(A(3,4)
     &   *A(4,1)-A(3,1)*A(4,4))+A(2,4)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))
      B(4,1) = d*(A(2,1)*(A(3,3)*A(4,2)-A(3,2)*A(4,3))+A(2,2)*(A(3,1)
     &   *A(4,3)-A(3,3)*A(4,1))+A(2,3)*(A(3,2)*A(4,1)-A(3,1)*A(4,2)))
      B(1,2) = d*(A(1,2)*(A(3,4)*A(4,3)-A(3,3)*A(4,4))+A(1,3)*(A(3,2)
     &   *A(4,4)-A(3,4)*A(4,2))+A(1,4)*(A(3,3)*A(4,2)-A(3,2)*A(4,3)))
      B(2,2) = d*(A(1,1)*(A(3,3)*A(4,4)-A(3,4)*A(4,3))+A(1,3)*(A(3,4)
     &   *A(4,1)-A(3,1)*A(4,4))+A(1,4)*(A(3,1)*A(4,3)-A(3,3)*A(4,1)))
      B(3,2) = d*(A(1,1)*(A(3,4)*A(4,2)-A(3,2)*A(4,4))+A(1,2)*(A(3,1)
     &   *A(4,4)-A(3,4)*A(4,1))+A(1,4)*(A(3,2)*A(4,1)-A(3,1)*A(4,2)))
      B(4,2) = d*(A(1,1)*(A(3,2)*A(4,3)-A(3,3)*A(4,2))+A(1,2)*(A(3,3)
     &   *A(4,1)-A(3,1)*A(4,3))+A(1,3)*(A(3,1)*A(4,2)-A(3,2)*A(4,1)))
      B(1,3) = d*(A(1,2)*(A(2,3)*A(4,4)-A(2,4)*A(4,3))+A(1,3)*(A(2,4)
     &   *A(4,2)-A(2,2)*A(4,4))+A(1,4)*(A(2,2)*A(4,3)-A(2,3)*A(4,2)))
      B(2,3) = d*(A(1,1)*(A(2,4)*A(4,3)-A(2,3)*A(4,4))+A(1,3)*(A(2,1)
     &   *A(4,4)-A(2,4)*A(4,1))+A(1,4)*(A(2,3)*A(4,1)-A(2,1)*A(4,3)))
      B(3,3) = d*(A(1,1)*(A(2,2)*A(4,4)-A(2,4)*A(4,2))+A(1,2)*(A(2,4)
     &   *A(4,1)-A(2,1)*A(4,4))+A(1,4)*(A(2,1)*A(4,2)-A(2,2)*A(4,1)))
      B(4,3) = d*(A(1,1)*(A(2,3)*A(4,2)-A(2,2)*A(4,3))+A(1,2)*(A(2,1)
     &   *A(4,3)-A(2,3)*A(4,1))+A(1,3)*(A(2,2)*A(4,1)-A(2,1)*A(4,2)))
      B(1,4) = d*(A(1,2)*(A(2,4)*A(3,3)-A(2,3)*A(3,4))+A(1,3)*(A(2,2)
     &   *A(3,4)-A(2,4)*A(3,2))+A(1,4)*(A(2,3)*A(3,2)-A(2,2)*A(3,3)))
      B(2,4) = d*(A(1,1)*(A(2,3)*A(3,4)-A(2,4)*A(3,3))+A(1,3)*(A(2,4)
     &   *A(3,1)-A(2,1)*A(3,4))+A(1,4)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)))
      B(3,4) = d*(A(1,1)*(A(2,4)*A(3,2)-A(2,2)*A(3,4))+A(1,2)*(A(2,1)
     &   *A(3,4)-A(2,4)*A(3,1))+A(1,4)*(A(2,2)*A(3,1)-A(2,1)*A(3,2)))
      B(4,4) = d*(A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))+A(1,2)*(A(2,3)
     &   *A(3,1)-A(2,1)*A(3,3))+A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)))
      RETURN 
      END 

合并一个单位阵求逆

      SUBROUTINE MATINV(A,B)
      DIMENSION C(4,8),A(4,4),B(4,4)
      DO 1  K = 1,4
      DO 2  J = 1,4
      C(K,J)  = A(K,J)
2     CONTINUE
      C(K,J+4)= 1.0
1     CONTINUE
      DO 10 K = 1,4
      DO 20 J = 8,k,-1
20    C(K,J)  = C(K,J)/C(K,K)
      DO 40 I = 1,5
      IF (I.NE.K) THEN
            do 30 j=10,k,-1
30          c(i,j)=c(i,j)-c(i,k)*c(k,j)
      ENDIF
40    CONTINUE
10    CONTINUE
      DO 3  K = 1,4
      DO 3  J = 1,4
3     B(K,J)  = C(K,J+4)
      RETURN
      END

MKL库函数

      SUBROUTINE MTINV(A,B)
      USE LAPACK95
      DIMENSION A(4,4),B(4,4)
      CALL GETRF(A,IP)
      CALL GETRI(A,IP)
      B = A
      RETURN
      END