Testing factorisation in Julia

FastBackslash.jl

using LinearAlgebra, DualNumbers, BenchmarkTools, SparseArrays, SuiteSparse

n = 1000 # size of J and y

# make y
a, b = rand(n), rand(n)
y = a + b * ε

# make J
spA = sprandn(n, n, 10/n) + I
i, j, vA = findnz(spA)
spB = sparse(i, j, randn(length(spA.nzval)))
i, j, vB = findnz(spB)
spJ = sparse(i, j, vA .+ ε * vB)
A, B = Matrix(spA), Matrix(spB)
J = A + B * ε

# make complex versions of J and spJ to compare as well
h = 1e-50
imJ = A + B * im * h
spimJ = spA + spB * im * h
cy = a + b * im * h

# Types that stores the required terms only
mutable struct DualFactors
   Af::LinearAlgebra.LU{Float64,Array{Float64,2}} # the factors of the real part
   B::Array{Float64,2} # the non-real part
end
mutable struct SparseDualFactors
   Af::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64} # the factors of the real part
   B::SparseMatrixCSC{Float64,Int64} # the non-real part
end

# factorization function acts as a constructor
fast_factorize(M::Array{Dual{Float64},2}) = DualFactors(factorize(realpart.(M)), dualpart.(M))
fast_factorize(M::SparseMatrixCSC{Dual{Float64},Int64}) = SparseDualFactors(factorize(realpart.(M)), dualpart.(M))

# overload backslash
function Base.:\(J::Union{DualFactors, SparseDualFactors}, y::Vector{Dual{Float64}})
    a, b = realpart.(y), dualpart.(y)
    Af, B = J.Af, J.B
    tmp = Af \ a
    return tmp + (Af \ (b - B * tmp)) * ε
end

# Check that the results are approximately the same
# 1. Factorize
Jf = factorize(J)
imJf = factorize(imJ)
spimJf = factorize(spimJ)
fastJf = fast_factorize(J)
fastspJf = fast_factorize(spJ)
# 2. Back-substitute
x = Jf \ y
imx = imJf \ cy
spimx = spimJf \ cy
fastx = fastJf \ y
fastspx = fastspJf \ y
# 3. Test for approximate equality of solution
realpart.(x) ≈ realpart.(fastx)
dualpart.(x) ≈ dualpart.(fastx)
realpart.(x) ≈ realpart.(fastspx)
dualpart.(x) ≈ dualpart.(fastspx)
realpart.(x) ≈ real.(imx)
dualpart.(x) ≈ imag.(imx / h)
realpart.(x) ≈ real.(spimx)
dualpart.(x) ≈ imag.(spimx / h)

# Start benchmark group
suite = BenchmarkGroup()
suite["fac"] = BenchmarkGroup(["factors", "factorize"])
suite["bs"] = BenchmarkGroup(["solve", "left-divide"])

# Benchmark group for the factorization
suite["fac"]["Dual"] = @benchmarkable Jf = factorize(J)
suite["fac"]["Complex"] = @benchmarkable imJf = factorize(imJ)
suite["fac"]["Sparse Complex"] = @benchmarkable spimJf = factorize(spimJ)
suite["fac"]["Fast Dual"] = @benchmarkable fastJf = fast_factorize(J)
suite["fac"]["Fast Sparse Dual"] = @benchmarkable fastspJf = fast_factorize(spJ)

# Benchmark group for the back-substitution
suite["bs"]["Dual"] = @benchmarkable x = Jf \ y
suite["bs"]["Complex"] = @benchmarkable imx = imJf \ cy
suite["bs"]["Sparse Complex"] = @benchmarkable spimx = spimJf \ cy
suite["bs"]["Fast Dual"] = @benchmarkable fastx = fastJf \ y
suite["bs"]["Fast Sparse Dual"] = @benchmarkable fastspx = fastspJf \ y

# Run the benchmarks
tune!(suite);
results = run(suite, verbose = true, seconds = 30)

# Show the benchmarks
results["fac"]
results["bs"]