实矩阵求逆的全选主元高斯-约当消去法

[Fortran] 纯文本查看 复制代码

PROGRAM test_14_Matrix_inv 
implicit none

!定义动态数组 

REAL,ALLOCATABLE :: X(:,:),L(:,:),U(:,:),y(:),z(:),inv_X(:,:) 

!定义矩阵的维数 

INTEGER N,I,J
write(*,*)'请输入方阵的维数N'
READ *,N 
ALLOCATE(X(N,N)) 
ALLOCATE(U(N,N)) 
ALLOCATE(L(N,N)) 
ALLOCATE(y(N)) 
ALLOCATE(z(N))
ALLOCATE(inv_X(N,N))
!输入矩阵X 

PRINT *,'请确保矩阵可逆!'
PRINT *,'请输入矩阵X',N,'行',N,'列' 

DO I=1,N 
PRINT *,'请输入第',I,'行' 
READ *,X(I,:) 
END DO 

!将矩阵X进行LU分解 

CALL Matrix_LU(X,N,L,U) 

!求解逆矩阵 

DO I=1,N 
CALL cal_y(L,y,N,I)  
CALL cal_z(U,y,z,N) 
DO J=1,N 
  inv_X(J,I)=z(J) 
END DO
END DO 

!打印结果 

DO I=1,N 
PRINT *,inv_X(I,:) 
END DO 
read(*,*)
END program test_14_Matrix_inv
!======================================================================! 
!        各子程序          ! 
!======================================================================! 
!子程序1=======LU分解 

SUBROUTINE Matrix_LU(X,N,L,U)
REAL L(N,N),U(N,N),X(N,N)
INTEGER I,N,J,K,M 
DO I=1,N
    U(1,I)=X(1,I) 
    L(I,1)=X(I,1)/U(1,1) 
END DO 

!求出U的第K行和U的第K列(K=2,3...N)

DO K=2,N
    DO J=K,N 
    U(K,J)=X(K,J)
    DO M=1,K-1    
        U(K,J)=U(K,J)-L(K,M)*U(M,J)
    END DO
END DO

IF (N==2) THEN 

  EXIT 

END IF

DO I=K+1,N      
  L(I,K)=X(I,K) 
  DO M=1,K-1
   L(I,K)=L(I,K)-L(I,M)*U(M,K) 
  ENDDO
  L(I,K)=L(I,K)/U(K,K) 
END DO 
END DO 

!L矩阵的主对角线全为1 

DO I=1,N 
L(I,I)=1
END DO 

!L为下三角,U为上三角

DO I=1,N     
DO J=1,N
  IF (I<J) THEN 
   L(I,J)=0 
  ENDIF 
END DO 
END DO  

DO I=1,N
    DO J=1,N 
  IF (I>J) THEN    
      U(I,J)=0 
  END IF 
END DO 
END DO

END SUBROUTINE Matrix_LU

!子程序2=======解Ly=b方程 

SUBROUTINE cal_y(L,y,N,K)
REAL y(N),L(N,N),B(N)
INTEGER I,J,N,K
DO I=1,N
B(I)=0 
END DO 
B(K)=1 
y(1)=B(1) 
DO I=2,N 
y(I)=B(I)
DO J=1,I-1 
  y(I)=y(I)-L(I,J)*y(J) 
END DO
END DO 
END SUBROUTINE cal_y 

!子程序3=======解Uz=y方程

SUBROUTINE cal_z(U,y,z,N)

REAL y(N),z(N),U(N,N) 
INTEGER I,J,N 
z(N)=y(N)/U(N,N) 
DO I=N-1,1  
    z(I)=y(I) 
DO J=I+1,N
  z(I)=z(I)-U(I,J)*z(J)
END DO 
z(I)=z(I)/U(I,I) 
END DO 
END SUBROUTINE cal_z 

[Fortran] 纯文本查看 复制代码

SUBROUTINE BRINV(A,N,L,IS,JS)
        DIMENSION A(N,N),IS(N),JS(N)
        DOUBLE PRECISION A,T,D
        L=1
        DO 100 K=1,N
          D=0.0
          DO 10 I=K,N
          DO 10 J=K,N
            IF (ABS(A(I,J)).GT.D) THEN
              D=ABS(A(I,J))
              IS(K)=I
              JS(K)=J
            END IF
10          CONTINUE
          IF (D+1.0.EQ.1.0) THEN
            L=0
            WRITE(*,20)
            RETURN
          END IF
20          FORMAT(1X,'ERR**NOT INV')
          DO 30 J=1,N
            T=A(K,J)
            A(K,J)=A(IS(K),J)
            A(IS(K),J)=T
30          CONTINUE
          DO 40 I=1,N
            T=A(I,K)
            A(I,K)=A(I,JS(K))
            A(I,JS(K))=T
40          CONTINUE
          A(K,K)=1/A(K,K)
          DO 50 J=1,N
            IF (J.NE.K) THEN
              A(K,J)=A(K,J)*A(K,K)
            END IF
50          CONTINUE
          DO 70 I=1,N
            IF (I.NE.K) THEN
              DO 60 J=1,N
                IF (J.NE.K) THEN
                  A(I,J)=A(I,J)-A(I,K)*A(K,J)
                END IF
60              CONTINUE
            END IF
70          CONTINUE
          DO 80 I=1,N
            IF (I.NE.K) THEN
              A(I,K)=-A(I,K)*A(K,K)
            END IF
80          CONTINUE
100        CONTINUE
        DO 130 K=N,1,-1
          DO 110 J=1,N
            T=A(K,J)
            A(K,J)=A(JS(K),J)
            A(JS(K),J)=T
110          CONTINUE
          DO 120 I=1,N
            T=A(I,K)
            A(I,K)=A(I,IS(K))
            A(I,IS(K))=T
120          CONTINUE
130        CONTINUE
        RETURN
        END

为啥不用OpenBLAS或者MKL呢？Windows下面直接用编译好的动态或者静态链接库，本站视频教程有解说
参考下面的代码，把B改成单位矩阵
http://www.nag.com/lapack-ex/node6.html