Naïve LU

module naivelu

using BenchmarkTools
using LinearAlgebra

"""
    naivelu!(a)

Naive LU decomposition implementation.

This function does not perform pivoting, that is no row or 
column exchanges are performed. 
"""
function naivelu!(a)
    n,  m = size(a);
    @assert n == m "Matrix must be square!"

    for col = 1:n-1 # for each column but the last
        piv = a[col, col]; # pivot
        col1=col+1;
        ipiv = 1.0 / piv
        @inbounds for r = col1:n # col of L    
            a[r, col] *= ipiv
        end
        # update submatrix A[col1:n,col1:n]
        @inbounds for c=col1:n
            nacolc = -a[col, c]
            @inbounds for r=col1:n
                a[r, c] += a[r, col]*nacolc; 
            end
        end
    end

    return a
end

function test()
    for i = 1:50
        N = 500
        a = rand(N, N);
        ap = deepcopy(a)
        F = lu!(ap, Val(false))
        a = naivelu!(a)
        U =  UpperTriangular(a)
        L =  LowerTriangular(a) - Diagonal(a) + UniformScaling(1.0)

        @assert  norm(F.L - L) < 1.0e-6
        @assert  norm(F.U - U) < 1.0e-6
    end
end

function timeit()

    N = 1000
    a = rand(N, N);
    @show N
    @btime lu!($a, Val(false))
    @btime naivelu!($a)

end

# using Profile  
# N = 5000                                                             
# a = rand(N, N);                                                         
# naivelu!(a);
# @profile naivelu!(a);                                                        
# Profile.print()                                                              
end

using Main.naivelu
# naivelu.test()
naivelu.timeit()

using PGFPlotsX
t1 = Table([:x => [1000, 2000, 3000, 4000, 6000, 8000, 10000], :y => [0.245, 2.094, 7.822, 19.846, 68.860, 165.281, 322.740]])
t2 = Table([:x => [1000, 2000, 3000, 4000, 6000, 8000, 10000], :y => [0.118, 1.161, 4.764, 14.134, 53.935, 133.058, 266.291]])
@pgf ax = LogLogAxis({
    height = "11cm",
    width = "14cm",
    xlabel = "Number of equations",
    ylabel = "Time [s]",
    grid="major",
    enlargelimits = false, 
    legend_pos  = "south east"
}, 
Plot({very_thick, mark = "o", color="green", enlargelimits = false}, t1), LegendEntry("Built-in"),
Plot({very_thick, mark = "x", color="black", enlargelimits = false}, t2), LegendEntry("Naïve")
)
display(ax)