Created
October 18, 2018 23:46
-
-
Save briochemc/594fa371676302290018ddcb8a60b341 to your computer and use it in GitHub Desktop.
Testing factorisation in Julia
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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"] | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment