Testing LU

testinglu.jl

module testinglu

using LinearAlgebra
import LinearAlgebra: lu!, generic_lufact!, BlasInt, checknonsingular

function better_generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot} = Val(true);
                         check::Bool = true) where {T,Pivot}
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    Pivot && (ipiv = Vector{BlasInt}(undef, minmn))
    @inbounds begin
        for k = 1:minmn
            # find index max
            kp = k
            if Pivot
                amax = abs(zero(T))
                for i = k:m
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            Pivot && (ipiv[k] = kp)
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:m
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = k+1:n
                nAkj = -A[k,j]
                for i = k+1:m
                    A[i,j] += A[i,k]*nAkj
                end
            end
        end
    end
    (!Pivot) && (ipiv = Vector{BlasInt}(1:minmn))
    check && checknonsingular(info)
    return LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
end

function timeit(func)
    @show Base.Threads.nthreads()
    Ns = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000]
    ts = Float64[]
    for N = Ns
        numtrials = (8000 / N)^2
        a = rand(N, N)
        b = deepcopy(a)
        t = 0.0
        for i = 1:numtrials
            copyto!(b, a)
            t = t + @elapsed  func(b, Val(true))
        end
        t = t / numtrials
        push!(ts, t)
    end
    return Ns, ts
end
                                                         
end

using Main.testinglu
using LinearAlgebra
using PGFPlotsX
# Obtained with my laptop:
# julia> versioninfo(verbose = true)                                                      
# Julia Version 1.0.1                                                                     
# Commit 0d713926f8 (2018-09-29 19:05 UTC)                                                
# Platform Info:                                                                          
#   OS: Windows (x86_64-w64-mingw32)                                                      
#       Microsoft Windows [Version 10.0.17134.345]                                        
#   CPU: Intel(R) Core(TM) i7-6650U CPU @ 2.20GHz:                                        
#               speed         user         nice          sys         idle          irq    
#        #1  2208 MHz   10071703            0     10006968     66909468      1159250  ticks                                                                                       
#        #2  2208 MHz    9976859            0      6879250     70131625       100687  ticks                                                                                       
#        #3  2208 MHz   12104953            0      8315828     66566953       112546  ticks                                                                                       
#        #4  2208 MHz   13647359            0      7742703     65597671        97421  ticks                                                                                       
                                                                                     
#   Memory: 15.927024841308594 GB (8624.875 MB free)                                      
#   Uptime: 167506.5131633 sec                                                            
#   Load Avg:  0.0  0.0  0.0                                                              
#   WORD_SIZE: 64                                                                         
#   LIBM: libopenlibm                                                                     
#   LLVM: libLLVM-6.0.0 (ORCJIT, skylake)                                              
# (Ns, ts) = testinglu.timeit(LinearAlgebra.lu!) = ([5, 10, 15, 20, 30, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000], [4.92825e-7, 1.3246e-6, 2.98111e-6, 7.30522e-5, 0.000232655, 0.000346878, 0.000783372, 0.00163609, 0.00296943, 0.00404206, 0.00481048, 0.00904308, 0.00774103, 0.0103748, 0.0130644, 0.0161573, 0.0939215])                 
# (Base.Threads).nthreads() = 1                                                           
# (Ns, ts) = testinglu.timeit(LinearAlgebra.generic_lufact!) = ([5, 10, 15, 20, 30, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000], [2.35853e-7, 7.96788e-7, 1.85203e-6, 4.16479e-6, 1.1344e-5, 4.22405e-5, 0.000313653, 0.00232984, 0.0076944, 0.0188079, 0.0345471, 0.071774, 0.0960256, 0.161923, 0.212274, 0.309943, 2.4843])                   
# (Base.Threads).nthreads() = 1                                                           
# (Ns, ts) = testinglu.timeit(testinglu.better_generic_lufact!) = ([5, 10, 15, 20, 30, 50, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000], [2.45576e-7, 6.3709e-7, 1.54438e-6, 3.0596e-6, 7.31294e-6, 2.26065e-5, 0.000147725, 0.00100786, 0.00300646, 0.00753821, 0.0137282, 0.026914, 0.0500318, 0.0755291, 0.114681, 0.142427, 1.45933])  
@show Ns, ts = testinglu.timeit(LinearAlgebra.lu!)
tblas = Table([:x => Ns, :y => ts])
@show Ns, ts = testinglu.timeit(LinearAlgebra.generic_lufact!)
tgeneric = Table([:x => Ns, :y => ts])
@show Ns, ts = testinglu.timeit(testinglu.better_generic_lufact!)
tbetter = Table([:x => Ns, :y => ts])

@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"}, tblas), LegendEntry("BLAS"),
Plot({very_thick, mark = "x", color="blue"}, tgeneric), LegendEntry("generic"),
Plot({very_thick, mark = "triangle", color="blue"}, tbetter), LegendEntry("better generic")
)
display(ax)