LU

subroutine LU(n, m, A, P, Y)
    integer, intent(in) :: n, m
    integer :: i, j, k, pivot, temp
    double precision :: pivotValue
    integer, dimension(:), intent(inout) :: P
    double precision, dimension(:), allocatable :: x
    double precision, dimension(:, :), intent(inout) :: A, Y
    allocate(x(n))
    do i = 1, n
        !searching for max
        pivotValue = 0.0
        pivot = -1
        do k = i, n
            if (abs(A(k, i)) > pivotValue) then
                pivotValue = abs(A(k, i))
                pivot = k
            end if
        end do
        !is matrix invertible?
        if (abs(pivotValue) < 1.0/10000) then
            print *, "matrix is invertiblle"
        end if
        !row swap
        !call dswap(1, P(i), 0, P(pivot), 0)
        temp = P(i)
        P(i) = P(pivot)
        P(pivot) = temp
        call dswap(n, A(i,1:n), 1, A(pivot, 1:n), 1)
        do j = i + 1, n
            A(j, i) = A(j, i) / A(i, i)
            do k = i + 1, n
                A(j, k) = A(j, k) - A(j, i) * A(i, k)
            end do
        end do
    end do
    !A = L + U - I
    !Gauss
    do k = 1, m
        !forward
        do i = 1, n
            x(i) = Y(P(i), k)
            x(i) = x(i) - ddot(i-1, A(i,1:i-1), 1, x(1:i-1), 1)
            !do j = 1, i - 1
            !    x(i) = x(i) - A(i, j) * x(j)
            !end do
        end do
        !backward
        do i = n, 1, -1
            if (i < n) then
                x(i) = x(i) - ddot(n-i, x(i+1:n), 1, A(i,i+1:n), 1)
            end if
            !do j = i + 1, n
            !    x(i) = x(i) - A(i, j) * x(j)
            !end do
            x(i) = x(i) / A(i, i)
        end do
        !answer
        do i = 1, n
            Y(i, k) = x(i)
        end do
    end do
    deallocate(x)
end subroutine

program lu_decomposition
    implicit none
    integer :: n, m, i, inp, outp
    integer, dimension(:), allocatable :: P, sP
    double precision, dimension(:,:), allocatable :: A, sA, tempY, Y, sY
    interface
        subroutine LU(n, m, A, P, Y)
            integer, intent(in) :: n, m
            integer :: i, j, k, pivot, temp
            double precision :: pivotValue
            integer, dimension(:), intent(inout) :: P
            double precision, dimension(:), allocatable :: x
            double precision, dimension(:, :), intent(inout) :: A, Y
        end subroutine LU
    end interface
    !input
    inp = 42
    outp = 43
    open(unit = inp, file = "input.txt")
    read(inp, *) n, m
    allocate(P(n))
    do i = 1, n
        P(i) = i
    end do
    allocate(sP(n))
    sP = P
    allocate(A(n, n))
    allocate(sA(n, n))
    allocate(Y(n, m))
    allocate(sY(n, m))
    allocate(tempY(m, n))
    read (inp, *) A
    A = transpose(A)
    sA = A
    read (inp, *) tempY
    Y = transpose(tempY)
    sY = Y
    close(inp)
    open(unit = outp, file = "output.txt")
    !LU decomposition
    call LU(n, m, A, P, Y)
    !LU output
    write(outp, *) "A"
    do i = 1, n
        write(outp, *) A(i, :)
    end do
    write(outp, *) "P"
    write(outp, *) P
    write(outp, *) "Y"
    do i = 1, n
        write(outp, "(1000F10.3)") Y(i, :)
    end do
    !sgesv
    call dgesv(n, m, sA, n, sP, sY, n, i)
    !sgesv output
    if (i .eq. 0) then
        print *, "successful sgesv execution"
    end if
    write(outp, *) "A by sgesv"
    do i = 1, n
        write(outp, *) sA(i, :)
    end do
    write(outp, *) "P by sgesv"
    write(outp, *) sP
    write(outp, *) "Y by sgesv"
    do i = 1, n
        write(outp, "(1000F10.3)") sY(i, :)
    end do
    close(outp)
    !memory
    deallocate(P)
    deallocate(sP)
    deallocate(A)
    deallocate(Y)
    deallocate(sY)
end program