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

[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