求逆矩阵的程序，结果只有最后一行是对的

little_kar · 发表于 2017-2-20 17:40:32

用下面这个程序求矩阵的逆，结果却总是只有最后一行是正确的，而前n-1行的数据永远全部是-4.3160208E+08。一下是我试验的矩阵、结果和程序：

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

001PROGRAM test_14_Matrix_inv 
002implicit none
003 
004!定义动态数组 
005 
006REAL,ALLOCATABLE :: X(:,:),L(:,:),U(:,:),y(:),z(:),inv_X(:,:)
007 
008!定义矩阵的维数 
009 
010INTEGER N,I,J
011write(*,*)'请输入方阵的维数N'
012READ *,N 
013ALLOCATE(X(N,N))
014ALLOCATE(U(N,N))
015ALLOCATE(L(N,N))
016ALLOCATE(y(N))
017ALLOCATE(z(N))
018ALLOCATE(inv_X(N,N))
019!输入矩阵X 
020 
021PRINT *,'请确保矩阵可逆!'
022PRINT *,'请输入矩阵X',N,'行',N,'列'
023 
024DO I=1,N 
025PRINT *,'请输入第',I,'行'
026READ *,X(I,:)
027END DO
028 
029!将矩阵X进行LU分解 
030 
031CALL Matrix_LU(X,N,L,U) 
032 
033!求解逆矩阵 
034 
035DO I=1,N 
036CALL cal_y(L,y,N,I)  
037CALL cal_z(U,y,z,N) 
038DO J=1,N 
039  inv_X(J,I)=z(J)
040END DO
041END DO
042 
043!打印结果 
044 
045DO I=1,N 
046PRINT *,inv_X(I,:)
047END DO
048read(*,*)
049END program test_14_Matrix_inv
050!======================================================================! 
051!        各子程序          ! 
052!======================================================================! 
053!子程序1=======LU分解 
054 
055SUBROUTINE Matrix_LU(X,N,L,U)
056REAL L(N,N),U(N,N),X(N,N)
057INTEGER I,N,J,K,M 
058DO I=1,N
059    U(1,I)=X(1,I)
060    L(I,1)=X(I,1)/U(1,1)
061END DO
062 
063!求出U的第K行和U的第K列(K=2,3...N)
064 
065DO K=2,N
066    DO J=K,N 
067    U(K,J)=X(K,J)
068    DO M=1,K-1   
069        U(K,J)=U(K,J)-L(K,M)*U(M,J)
070    END DO
071END DO
072 
073IF (N==2) THEN
074 
075  EXIT
076 
077END IF
078 
079DO I=K+1,N      
080  L(I,K)=X(I,K)
081  DO M=1,K-1
082   L(I,K)=L(I,K)-L(I,M)*U(M,K)
083  ENDDO
084  L(I,K)=L(I,K)/U(K,K)
085END DO
086END DO
087 
088!L矩阵的主对角线全为1 
089 
090DO I=1,N 
091L(I,I)=1
092END DO
093 
094!L为下三角,U为上三角
095 
096DO I=1,N     
097DO J=1,N
098  IF (I<J) THEN
099   L(I,J)=0
100  ENDIF
101END DO
102END DO 
103 
104DO I=1,N
105    DO J=1,N 
106  IF (I>J) THEN   
107      U(I,J)=0
108  END IF
109END DO
110END DO
111 
112END SUBROUTINE Matrix_LU
113 
114!子程序2=======解Ly=b方程 
115 
116SUBROUTINE cal_y(L,y,N,K)
117REAL y(N),L(N,N),B(N)
118INTEGER I,J,N,K
119DO I=1,N
120B(I)=0
121END DO
122B(K)=1
123y(1)=B(1)
124DO I=2,N 
125y(I)=B(I)
126DO J=1,I-1
127  y(I)=y(I)-L(I,J)*y(J)
128END DO
129END DO
130END SUBROUTINE cal_y 
131 
132!子程序3=======解Uz=y方程
133 
134SUBROUTINE cal_z(U,y,z,N)
135 
136REAL y(N),z(N),U(N,N)
137INTEGER I,J,N 
138z(N)=y(N)/U(N,N)
139DO I=N-1,1 
140    z(I)=y(I)
141DO J=I+1,N
142  z(I)=z(I)-U(I,J)*z(J)
143END DO
144z(I)=z(I)/U(I,I)
145END DO
146END SUBROUTINE cal_z

pasuka · 发表于 2017-2-20 20:57:42

本帖最后由 pasuka 于 2017-2-20 20:59 编辑

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

fcode · 发表于 2017-2-21 22:49:09

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

01SUBROUTINE cal_z(U,y,z,N)
02 
03REAL y(N),z(N),U(N,N)
04INTEGER I,J,N 
05z(N)=y(N)/U(N,N)
06DO I=N-1,1,-1 !//此处步长为 -1
07    z(I)=y(I)
08DO J=I+1,N
09  z(I)=z(I)-U(I,J)*z(J)
10END DO
11z(I)=z(I)/U(I,I)
12END DO
13END SUBROUTINE cal_z

V_Sc · 发表于 2017-3-1 10:59:55

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

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

01SUBROUTINE BRINV(A,N,L,IS,JS)
02        DIMENSION A(N,N),IS(N),JS(N)
03        DOUBLE PRECISION A,T,D
04        L=1
05        DO 100 K=1,N
06          D=0.0
07          DO 10 I=K,N
08          DO 10 J=K,N
09            IF (ABS(A(I,J)).GT.D) THEN
10              D=ABS(A(I,J))
11              IS(K)=I
12              JS(K)=J
13            END IF
1410          CONTINUE
15          IF (D+1.0.EQ.1.0) THEN
16            L=0
17            WRITE(*,20)
18            RETURN
19          END IF
2020          FORMAT(1X,'ERR**NOT INV')
21          DO 30 J=1,N
22            T=A(K,J)
23            A(K,J)=A(IS(K),J)
24            A(IS(K),J)=T
2530          CONTINUE
26          DO 40 I=1,N
27            T=A(I,K)
28            A(I,K)=A(I,JS(K))
29            A(I,JS(K))=T
3040          CONTINUE
31          A(K,K)=1/A(K,K)
32          DO 50 J=1,N
33            IF (J.NE.K) THEN
34              A(K,J)=A(K,J)*A(K,K)
35            END IF
3650          CONTINUE
37          DO 70 I=1,N
38            IF (I.NE.K) THEN
39              DO 60 J=1,N
40                IF (J.NE.K) THEN
41                  A(I,J)=A(I,J)-A(I,K)*A(K,J)
42                END IF
4360              CONTINUE
44            END IF
4570          CONTINUE
46          DO 80 I=1,N
47            IF (I.NE.K) THEN
48              A(I,K)=-A(I,K)*A(K,K)
49            END IF
5080          CONTINUE
51100        CONTINUE
52        DO 130 K=N,1,-1
53          DO 110 J=1,N
54            T=A(K,J)
55            A(K,J)=A(JS(K),J)
56            A(JS(K),J)=T
57110          CONTINUE
58          DO 120 I=1,N
59            T=A(I,K)
60            A(I,K)=A(I,IS(K))
61            A(I,IS(K))=T
62120          CONTINUE
63130        CONTINUE
64        RETURN
65        END

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

回帖奖励 +20

评分

水王勋章

元老勋章

热心勋章

美女勋章

星光勋章

新人勋章

贡献勋章

管理勋章

帅哥勋章

爱心勋章

规矩勋章