[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