# L-10 MCS 471 Wed 14 Sep 2022 : plufac.jl
# Illustrates the formulas for LU factorization with row pivoting.

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.3e", f)

using LinearAlgebra

"""
    lufac(mat::Array{Float64,2}, verbose::Bool=true) 

returns the matrices L, U, and P in an LU factorization of mat,
with row pivoting.  If verbose, writes in every step the matrix.

Example:

    A = rand(4,4)
    L, U, P = lufac(A)
"""
function lufac(mat::Array{Float64,2}, verbose::Bool=true) 
    nbrows, nbcols = size(mat)
    low = Matrix{Float64}(I,nbrows, nbcols)
    upp = zeros(nbrows, nbcols)
    piv = [i for i=1:nbrows]
    wrk = deepcopy(mat)
    for j=1:nbcols
        valmax, idxmax = abs(wrk[j,j]), j
        for i=j+1:nbrows
            if abs(wrk[i,j]) > valmax
                valmax, idxmax = abs(wrk[i,j]), i
            end
        end
        if verbose
            print("Step ", j)
            println(", the pivot row : ", idxmax)
        end
        if idxmax != j
            wrk[idxmax,:], wrk[j,:] = wrk[j,:], wrk[idxmax,:]
            piv[idxmax], piv[j] = piv[j], piv[idxmax]
            if verbose
                println("The matrix after swap :")
                show(stdout, "text/plain", wrk); println("")
            end
        end
        for i=j+1:nbrows
            wrk[i, j] = wrk[i,j]/wrk[j,j]
            for k=j+1:nbcols
                wrk[i, k] = wrk[i, k] - wrk[i, j]*wrk[j, k]
            end
        end
        if verbose
            println("The matrix after the update :")
            show(stdout, "text/plain", wrk); println("")
        end           
    end
    for i=1:nbrows
        for j=1:i-1
            low[i, j] = wrk[i, j]
        end
        for j=i:nbcols
            upp[i, j] = wrk[i, j]
        end
    end
    return [low, upp, piv]
end

"""
Prompts for a dimension and generates a random matrix,
returned by the function.
"""
function random_example()
    print("Give the dimension : ")
    line = readline(stdin)
    dim = parse(Int, line)
    mat = rand(dim, dim)
    return mat
end

"""
Returns the example used in lecture 10.
"""
function lecture_example()
    mat = [ 9.229e-2, -1.324e+0,  1.976e+0,
           -6.501e-1,  1.201e+0, -3.308e-1,
            2.245e+0, -1.265e+0, -1.277e+0]
    mat = reshape(mat, (3, 3))
    mat = permutedims(mat)
    return mat
end

"""
Prompts the user for a dimension and then
generates a random matrix to test the formulas
for a LU factorization.
"""
function main()
    print("Random matrix ? (y/n) ")
    line = readline(stdin)
    if line[1] == 'y'
        mat = random_example()
    else
        mat = lecture_example()
    end
    println("A matrix :")
    show(stdout, "text/plain", mat); println("");
    low, upp, piv = lufac(mat)
    println("The output of the LU factorization :")
    println("The pivots : ", piv)
    println("The lower triangular matrix :")
    show(stdout, "text/plain", low); println("");
    println("The upper triangular matrix :")
    show(stdout, "text/plain", upp); println("");
    println("Verification of the output :")
    prod = low*upp
    println("The product of lower with upper :")
    show(stdout, "text/plain", prod); println("");
    println("The permuted input matrix :")
    show(stdout, "text/plain", mat[piv, :]); println("");
    difference = mat[piv, :] - prod
    println("The difference : ", norm(difference))
end

main()