Author:RainMan
!2009-12-20
!LU分解法解线性方程组AX=B
!该程序从文件输入，格式为
!矩阵的维数是 3
!A矩阵
!4 2 1
!8 7 4
!16 17 15
!B矩阵
!11 34 95
PROGRAM solve_equation
    REAL,DIMENSION(:,:),ALLOCATABLE::A
    REAL,DIMENSION(:),ALLOCATABLE::B
    CHARACTER*(80):: inputfile = 'dengjin.txt'
    INTEGER DIM
    CALL GET_DIM(inputfile,DIM)
    ALLOCATE(A(DIM,DIM),B(DIM))
    CALL INPUT_MATRIX(A,B,DIM,inputfile)
    CALL LU(A,DIM)
    CALL solve(A,B,B,DIM)
    PRINT*,B
END

SUBROUTINE GET_DIM(inputfile,n)
    CHARACTER*(*) inputfile
    CHARACTER*(80) str
    INTEGER n
    OPEN(10,FILE = inputfile)
    READ(10,*)str,n
    CLOSE(10)
END

SUBROUTINE INPUT_MATRIX(A,B,DIM,inputfile)
    CHARACTER*(*) inputfile
    INTEGER DIM
    INTEGER A(DIM,DIM),B(DIM)
   
    OPEN(20,FILE=inputfile)
    READ(20,*)
    READ(20,*)
    READ(20,*) ((A(I,J),J=1,DIM),I=1,DIM)
    READ(20,*)
    READ(20,*) (B(I),I=1,DIM)
    CLOSE(20)
END

SUBROUTINE LU(A,DIM)
    INTEGER DIM
    REAL A(DIM,DIM)
   
    DO I=2,DIM
        DO J=I,DIM
            A(J,I-1) = A(J,I-1)*1.0/A(I-1,I-1)
            m = A(J,I-1)
            DO K=I,DIM
                A(J,K) = A(J,K)-m*A(I-1,K)
            ENDDO
        ENDDO
    ENDDO
END

SUBROUTINE solve(A,B,X,DIM)
    INTEGER:: DIM
    REAL A(DIM,DIM),B(DIM),X(DIM)
!----------LY = B
    DO I=2,DIM
        DO J=1,I-1
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
    ENDDO

    B(DIM) = B(DIM)/A(DIM,DIM)
!---------UX=Y Y stores in B,same as X
    DO I=DIM-1,1,-1
        DO J=I+1,DIM
            B(I) = B(I) - A(I,J)*B(J)
        ENDDO
        B(I) = B(I)/A(I,I)
    ENDDO
END