Solve Tridiagonal Linear Systems using CROUT Factorization using FORTRAN 90/95 [LU FACTORIZATION FOR TRIDIAGONAL SYSTEM]

A02.F90

!Name      : LU FACTORIZATION FOR TRIDIAGONAL SYSTEM
!Author    : KAWSER
!Blog      : http://blog.kawser.org
!Created   : 22/05/2014 7:48 AM
!
!Short URL : http://goo.gl/Uu9dNV
!
!Purpose   : We'll solve the system using CROUT Factorization
!            for Tridiagonal Linear Systems
!
!
!Algorithm Reference : NUMERIACAL ANALYSIS
!                        By RICHARD L. BURDEN
!                           J.DOUGLAS FAIRES
!                      Page No - 408
!
!Sample input file "A02.txt"
!4
!2 -1 0 0 1
!-1 2 -1 0 0
!0 -1 2 -1 0
!0 0 -1 2 1







PROGRAM A02
    IMPLICIT NONE

    INTEGER::N !DIMENSION OF THE MATRIX

    REAL,ALLOCATABLE,DIMENSION(:,:)::A !AUGMENTED MATRIX
    REAL,ALLOCATABLE,DIMENSION(:,:)::L !LOWER TRIANGULAR MATRIX
    REAL,ALLOCATABLE,DIMENSION(:,:)::U !UPPER TRIANGULAR MATRIX
    REAL,ALLOCATABLE,DIMENSION(:)::X !SOLUTION OF THE SYSTEM

    LOGICAL::SOLVE_TRIDIAGONAL_SYSTEM , R
    INTEGER::ROW,COLUMN


    OPEN(10 , FILE = "A02.txt") !OPENING INPUT FILE
    OPEN(20 , FILE = "A02_OUT.txt") !OPENING OUTPUT FILE

    READ(10,*) N !READING THE DIMENSION OF THE MATRIX FROM THE FILE

    ALLOCATE( A(N,N+1) , L(N,N+1) , U(N,N+1) , X(N) ) !ALLOCATING MEMORY

    !READING THE AUGMENTED MATRIX
    READ(10,*) (  ( A(ROW,COLUMN) , COLUMN=1 , N+1 ) , ROW=1 , N)
    100 FORMAT(1X,F10.3) !THIS FORMATTING STYLE IS USED TO FORMAT THE MATRIX PROPERLY

    WRITE(20,*) "#-------------------- AUGMENTED MATRIX A|b ----------------------------------------#"
    DO ROW = 1,N
        DO COLUMN = 1,N+1
            WRITE(20,100,ADVANCE='NO') A(ROW,COLUMN)
        END DO
        WRITE(20,*) !MOVE WRITE CURSOR TO NEW LINE
    END DO

    !CALLING TRIDIAGONAL SYSTEM SOLVER
    R = SOLVE_TRIDIAGONAL_SYSTEM(N,A,L,U,X)
    IF( R .EQV. .TRUE. ) THEN
        WRITE(20,*) "LOWER TRIANGULA MATRIX"
        DO ROW = 1,N
            DO COLUMN = 1,N
                WRITE(20,100,ADVANCE='NO') L(ROW,COLUMN)
            END DO
            WRITE(20,*) !MOVE WRITE CURSOR TO NEW LINE
        END DO

        WRITE(20,*) "UPPER TRIANGULA MATRIX"
        DO ROW = 1,N
            DO COLUMN = 1,N
                WRITE(20,100,ADVANCE='NO') U(ROW,COLUMN)
            END DO
            WRITE(20,*) !MOVE WRITE CURSOR TO NEW LINE
        END DO

        WRITE(20,*) "SOLUTION X"
        DO ROW=1,N
            WRITE(20,100,ADVANCE="NO") X(ROW)
        END DO
        WRITE(20,*) !MOVE WRITE CURSOR TO NEW LINE
    ELSE
        WRITE(20,*) "SORRY THE SYSTEM CAN NOT BE SOLVED USING CROUT FACTORIZATION"
    END IF

    DEALLOCATE( A , L , U , X ) !DEALLOCATING ALLOCATED MEMORY
    CLOSE(10) !CLOSING INPUT FILE
    CLOSE(20) !CLOSING OUTPUT FILE
END PROGRAM





LOGICAL FUNCTION SOLVE_TRIDIAGONAL_SYSTEM( N , A , L , U ,X )
    IMPLICIT NONE

    INTEGER,INTENT(IN)::N

    REAL,INTENT(IN),DIMENSION(N,N+1)::A !AUGMENTED MATRIX
    REAL,INTENT(OUT),DIMENSION(N,N+1)::L !LOWER TRIANGUAL MATRIX
    REAL,INTENT(OUT),DIMENSION(N,N+1)::U !UPPER TRIANGUAL MATRIX
    REAL,INTENT(OUT),DIMENSION(N)::X !SOLUTION OF THE SYSTEM

    REAL,DIMENSION(N)::Z
    INTEGER::I !LOOP VARIABLES

    !SOLVE Lz = b
    L(1,1) = A(1,1)
    IF ( L(1,1) == 0 ) THEN
        !CAN NOT BE SOLVED USING CROUT FACTORIZATION
        SOLVE_TRIDIAGONAL_SYSTEM = .FALSE.
        RETURN
    END IF
    U(1,1) = 1.0
    U(1,2) = A(1,2) / L(1,1)
    Z(1) = A(1,N+1) / L(1,1)

    DO I=2,N-1
        L(I,I-1) = A(I,I-1)
        L(I,I) = A(I,I) - L(I,I-1) * U(I-1,I)

        IF ( L(I,I) == 0 ) THEN
            !CAN NOT BE SOLVED USING CROUT FACTORIZATION
            SOLVE_TRIDIAGONAL_SYSTEM = .FALSE.
            RETURN
        END IF

        U(I,I) = 1.0
        U(I,I+1) = A(I,I+1) / L(I,I)
        Z(I) = ( A(I,N+1) - L(I,I-1) * Z(I-1) ) / L(I,I)
    END DO

    U(N,N) = 1.0
    L(N,N-1) = A(N,N-1)
    L(N,N) = A(N,N) - L(N,N-1) * U(N-1,N)
    IF ( L(N,N) == 0 ) THEN
        !CAN NOT BE SOLVED USING CROUT FACTORIZATION
        SOLVE_TRIDIAGONAL_SYSTEM = .FALSE.
        RETURN
    END IF

    Z(N) = ( A(N,N+1) - L(N,N-1) * Z(N-1) ) / L(N,N)

    !SOLVE Ux = z
    X(N) = Z(N)
    DO I = N-1,1,-1
        X(I) = Z(I) - U(I,I+1) * X(I+1)
    END DO

    SOLVE_TRIDIAGONAL_SYSTEM = .TRUE.
    RETURN
END FUNCTION