using BenchmarkTools, Random, VectorizationBase, Statistics
using LinearAlgebra, SparseArrays, LinearSolve, Sparspak
import Pardiso
using Plots
using MatrixDepot
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
# Why do I need to set this ?
BenchmarkTools.DEFAULT_PARAMETERS.samples = 10
algs = [
UMFPACKFactorization(),
KLUFactorization(),
MKLPardisoFactorize(),
SparspakFactorization(),
]
algnames = ["UMFPACK", "KLU", "Pardiso", "Sparspak"]
algnames_transpose = reshape(algnames, 1, length(algnames))
cols = [:red, :blue, :green, :magenta, :turqoise] # one color per alg
# matrices = ["HB/1138_bus", "HB/494_bus", "HB/662_bus", "HB/685_bus", "HB/bcsstk01", "HB/bcsstk02", "HB/bcsstk03", "HB/bcsstk04", "HB/bcsstk05", "HB/bcsstk06", "HB/bcsstk07", "HB/bcsstk08", "HB/bcsstk09", "HB/bcsstk10", "HB/bcsstk11", "HB/bcsstk12", "HB/bcsstk13", "HB/bcsstk14", "HB/bcsstk15", "HB/bcsstk16"]
allmatrices_md = listnames("*/*")
@info "Total number of matrices: $(allmatrices_md.content[1].rows)"
times = fill(NaN, length(allmatrices_md.content[1].rows), length(algs))
percentage_sparsity = fill(NaN, length(allmatrices_md.content[1].rows))
spaced_out_sparsity = fill(NaN, length(allmatrices_md.content[1].rows))
matrix_size = fill(NaN, length(allmatrices_md.content[1].rows))
bandedness_five = fill(NaN, length(allmatrices_md.content[1].rows))
bandedness_ten = fill(NaN, length(allmatrices_md.content[1].rows))
bandedness_twenty = fill(NaN, length(allmatrices_md.content[1].rows))
function compute_bandedness(A, bandwidth)
n = size(A, 1)
total_band_positions = 0
non_zero_in_band = 0
bandwidth = bandwidth
for r in 1:n
for c in 1:n
if abs(r - c) <= bandwidth
total_band_positions += 1 # This position belongs to the band
if A[r, c] != 0
non_zero_in_band += 1 # This element is non-zero in the band
end
end
end
end
percentage_filled = non_zero_in_band / total_band_positions * 100
return percentage_filled
end
compute_bandedness (generic function with 1 method)
for z in 1:length(allmatrices_md.content[1].rows)
try
matrix = allmatrices_md.content[1].rows[z]
matrix = string(matrix[1])
currMTX = matrix
rng = MersenneTwister(123)
A = mdopen(currMTX).A
A = convert(SparseMatrixCSC, A)
n = size(A, 1)
mtx_copy = copy(A)
@info "$n × $n"
n > 100 && error("Skipping too large matrices")
## COMPUTING SPACED OUT SPARSITY
rows, cols = size(mtx_copy)
new_rows = div(rows, 2)
new_cols = div(cols, 2)
condensed = zeros(Int, new_rows, new_cols)
while size(mtx_copy, 1) > 32 || size(mtx_copy, 2) > 32
rows, cols = size(mtx_copy)
new_rows = div(rows, 2)
new_cols = div(cols, 2)
condensed = sparse(zeros(Int, new_rows, new_cols))
for r in 1:2:rows-1
for c in 1:2:cols-1
block = mtx_copy[r:min(r+1, rows), c:min(c+1, cols)]
condensed[div(r-1, 2) + 1, div(c-1, 2) + 1] = (length(nonzeros(block)) >= 2) ? 1 : 0
end
end
mtx_copy = condensed
end
## COMPUTING FACTORIZATION TIME
b = rand(rng, n)
u0 = rand(rng, n)
for j in 1:length(algs)
bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem(copy($A),
copy($b);
u0 = copy($u0),
alias_A = true,
alias_b = true))
times[z,j] = bt
end
bandedness_five[z] = compute_bandedness(A, 5)
bandedness_ten[z] = compute_bandedness(A, 10)
bandedness_twenty[z] = compute_bandedness(A, 20)
percentage_sparsity[z] = length(nonzeros(A)) / n^2
spaced_out_sparsity[z] = length(nonzeros(mtx_copy)) * percentage_sparsity[z]
matrix_size[z] = n
#=
p = bar(algnames, times[z, :];
ylabel = "Time/s",
yscale = :log10,
title = "Time on $(currMTX)",
fmt = :png,
legend = :outertopright)
display(p)
=#
println("successfully factorized $(currMTX)")
catch e
matrix = allmatrices_md.content[1].rows[z]
matrix = string(matrix[1])
currMTX = matrix
println("$(currMTX) failed to factorize.")
println(e)
end
end
percentage_sparsity = percentage_sparsity[.!isnan.(percentage_sparsity)]
spaced_out_sparsity = spaced_out_sparsity[.!isnan.(spaced_out_sparsity)]
spaced_out_sparsity = replace(spaced_out_sparsity, 0 => 1e-10)
bandedness_five = bandedness_five[.!isnan.(bandedness_five)]
bandedness_five = replace(bandedness_five, 0 => 1e-10)
bandedness_ten = bandedness_ten[.!isnan.(bandedness_ten)]
bandedness_ten = replace(bandedness_ten, 0 => 1e-10)
bandedness_twenty = bandedness_twenty[.!isnan.(bandedness_twenty)]
bandedness_twenty = replace(bandedness_twenty, 0 => 1e-10)
matrix_size = matrix_size[.!isnan.(matrix_size)]
nanrows = any(isnan, times; dims=2)
times = times[.!vec(nanrows), :]