Fortran Coder

查看: 5908|回复: 3
打印 上一主题 下一主题

[线性代数] 求逆矩阵的程序,结果只有最后一行是对的

[复制链接]

35

帖子

11

主题

0

精华

熟手

F 币
129 元
贡献
118 点
跳转到指定楼层
楼主
发表于 2017-2-20 17:40:32 | 显示全部楼层 |只看大图 回帖奖励 |倒序浏览 |阅读模式
用下面这个程序求矩阵的逆,结果却总是只有最后一行是正确的,而前n-1行的数据永远全部是-4.3160208E+08。一下是我试验的矩阵、结果和程序:
[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 


1487583424(1).jpg (1.89 KB, 下载次数: 328)

要求的矩阵

要求的矩阵

1487583392(1).jpg (33.02 KB, 下载次数: 330)

得出的结果

得出的结果
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-1 06:58

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表