ReadMe.md

Backsolve after sparse LU: MATLAB vs Julia

See JuliaLang/julia issue #31105.

benchmark_Julia.jl

using LinearAlgebra, SparseArrays, SuiteSparse, BenchmarkTools, MAT
M = matread("test_matrix.mat")["M"] ;
x = matread("test_matrix.mat")["x"] ;
Mf = factorize(M) ;
@btime Mf = factorize(M) ;
@btime M \ x ;
@btime Mf \ x ;
@btime Mf \ [x 2x 3x 4x] ;

benchmark_MATLAB.m

n = 20000;
M = sprand(n, n, 20/n) + speye(n);
x = rand(n, 1);
tic; Mf = linfactor(M); toc;
tic; M \ x; toc;
tic; linfactor(Mf, x); toc;
tic; linfactor(Mf, [x 2*x 3*x 4*x]); toc;
save test_matrix.mat M x
# Error check
sol = linfactor(Mf, x) ;
norm(x)
norm(x - M * sol) / norm(x)
norm(x - M * sol)

Error_check.jl

# Error check
using LinearAlgebra, SparseArrays, SuiteSparse, BenchmarkTools, MAT
M = matread("test_matrix.mat")["M"] ;
x = matread("test_matrix.mat")["x"] ;
Mf = factorize(M) ;
sol1 = Mf \ x ;
sol2 = M \ x ;
norm(x)
norm(x - M * sol1) / norm(x)
norm(x - M * sol1)
norm(x - M * sol2) / norm(x)
norm(x - M * sol2)

Error_check.m

load test_matrix.mat M x
Mf = linfactor(M) ;
sol1 = linfactor(Mf, x) ;
sol2 = M \ x ;
norm(x)
norm(x - M * sol1) / norm(x)
norm(x - M * sol1)
norm(x - M * sol2) / norm(x)
norm(x - M * sol2)