PROGRAM ca5
    IMPLICIT NONE
    INTEGER :: i, j, k, n
    REAL :: a(5,5), b(5,1), l(5,5), u(5,5), x(5,1), y(5,1)

    !declaration
    n = 5
    a = RESHAPE( (/ 14.0, 14.0, -9.0, 3.0, 5.0, &
          & 14.0, 52.0, -15.0, 2.0, -32.0, &
          & -9.0, -15.0, 36.0, -5.0, 16.0, &
          & 3.0, 2.0, -5.0, 47.0, 49.0, &
          & -5.0, -32.0, 16.0, 49.0, 79.0 /), (/ 5, 5 /) )
    b = RESHAPE( (/ -15.0, -100.0, 106.0, 329.0, 463.0 /), (/ 5, 1 /) )
    DO i = 1, n
        DO j = 1, n
            l(i,j) = 0
            u(i,j) = 0
            x(i,j) = 0
            y(i,j) = 0
        ENDDO
        l(i,i) = 1
    ENDDO

    !LU, alternately solve row for upper and column for lower
    DO i = 1, n
        DO j = i, n
            u(i,j) = a(i,j)
            DO k = 1, i - 1
                u(i,j) = u(i,j) - u(k,j) * l(i,k)
            ENDDO
        ENDDO

        DO j = i + 1, n
            l(j,i) = a(j,i)
            DO k = 1, i - 1
                l(j,i) = l(j,i) - l(j,k) * u(k,i)
            ENDDO
            l(j,i) = l(j,i) / u(i,i)
        ENDDO
    ENDDO

    !find y
    DO i = 1, n
        y(i,1) = b(i,1)
        DO j = 1, i - 1
            y(i,1) = y(i,1) - l(i,j) * y(j,1)
        ENDDO
    ENDDO

    !find x
    DO i = 0, n - 1
        x(n-i,1) = y(n-i,1)
        DO j = 1, i
            x(n-i,1) = x(n-i,1) - u(n-i,n-j+1) * x(n-j+1,1)
        ENDDO
        x(n-i,1) = x(n-i,1) / u(n-i,n-i)
    ENDDO

    CALL print_mat('A',a)
    CALL print_mat('B',b)
    CALL print_mat('L',l)
    CALL print_mat('U',u)
    CALL print_mat('y',y)
    CALL print_mat('x',x, .TRUE.)
CONTAINS
    SUBROUTINE print_mat(name, mat, rounding)
        CHARACTER(1) :: name
        REAL, INTENT(IN) :: mat(:,:)
        LOGICAL, OPTIONAL :: rounding

        DO i = 1, 85
            WRITE(*, '(A)', ADVANCE = "NO") '+'
        ENDDO
        WRITE(*,*) name
        DO i = 1, n
            DO j = 1, INT(SIZE(mat)/n)
                WRITE(*,"(F10.5)", ADVANCE = "NO") mat(i,j)
                IF(PRESENT(rounding).AND.(rounding)) THEN
                    WRITE(*,*) "=~", NINT(mat(i,j))
                ENDIF
            ENDDO
            WRITE(*,*) ""
        ENDDO
    ENDSUBROUTINE
ENDPROGRAM ca5