LU factorization

program LU
    implicit none
    integer, parameter:: n=3
    integer:: i, j, k
    real:: A(n,n), b(n), x(n), y(n), L(n,n), U(n,n), s

    !reading
    open(11, file="in.txt")

    do i=1,n
        read(11,*) (A(i,j), j=1,n)
    end do

    read(11,*) (b(i), i=1,n)

    close(11)

    ! initializing 1st row of U and 1st column of L
    L(1,1)=1.0
    U(1,1)=A(1,1)

    if(A(1,1)==0.0) then
        write(*,*) "Impossible to find solution"
        stop
    end if

    do i=2,n
        U(1,i)=A(1,i)/L(1,1)
        L(i,1)=A(i,1)/U(1,1)
    end do


    !
    do i=2,n-1

        L(i,i)=1.0
        s=0.0
        do k=1,i-1
            s=s+L(i,k)*U(k,i)
        end do
        U(i,i)=A(i,i)-s

        if(A(i,i)==0.0) then
            write(*,*) "Impossible to find solution"
            stop
        end if

        do j=i+1,n

            s=0.0
            do k=1,i-1
                s=s+L(i,k)*U(k,j)
            end do
            U(i,j)=A(i,j)-s

            s=0.0
            do k=1,i-1
                s=s+L(j,k)*U(k,i)
            end do
            L(j,i)=(A(j,i)-s)/U(i,i)

        end do

    end do

    L(n,n)=1.0
    s=0.0
    do k=1,n-1
        s=s+L(n,k)*U(k,n)
    end do
    U(n,n)=A(n,n)-s


    open(12,file="out.txt")

    write(12,*)"Matrix L"
    do i=1,n
        write(12,*) (L(i,j), j=1,n)
    end do

    write(12,*)""
    write(12,*)"Matrix U"
    do i=1,n
        write(12,*) (U(i,j), j=1,n)
    end do

    ! forward substitution
    y(1)=b(1)
    do i=2,n
        s=0.0
        do j=1,i-1
            s=s+L(i,j)*y(j)
        end do
        y(i)=b(i)-s
    end do


    ! backward substitution
    x(n)=y(n)/U(n,n)
    do i=n-1,1,-1
        s=0.0
        do j=n,i+1,-1
            s=s+U(i,j)*x(j)
        end do
        x(i)=(y(i)-s)/U(i,i)
    end do

    write(12,*)""
    write(12,*)"solution is: "
    write(12,*) (x(i), i=1,n)

    close(12)

end program LU