ANN: Sparspak.jl 0.3.0 works with ForwardDiff.Dual and other number types

Package Announcements
1

This is a follow-up to the Announcement of Sparspak.jl, a sparse direct solver package resulting from the translation of the Sparspak Fortran90 library to Julia performed by @PetrKryslUCSD. With consent of the Fortran code authors Alan George and Joseph Liu, Sparspak.jl is released under an MIT license.

Starting with version 0.3.0, the internal calls to BLAS/LAPACK for number types not supported by the corresponding binary libraries are dispatched to specifically tailored, non-allocating implementations in Julia. This allows to auto-differentiate through a sparse matrix solve.

Here is an example:

using LinearAlgebra
using Sparspak
using Sparspak.SparseSolver: SparseSolver, solve!
using Sparspak.SpkProblem: Problem, insparse!, infullrhs!
using SparseArrays
using ForwardDiff


function tridiagonal(p,n)
    T=typeof(p)
    a=fill(T(-0.1),n-1)
    b=1.0.+p.*collect(1:n)
    Tridiagonal(a,b,a)
end
    

f(x)=sum(tridiagonal(x,20)\ones(20))
df(x)=ForwardDiff.derivative(f,x)

function g(x)
    n=20
    M=sparse(tridiagonal(x,n))
    pr = Problem(n,n,nnz(M),zero(x))
    insparse!(pr, M)
    infullrhs!(pr,ones(n))
    solve!(SparseSolver(pr))    
    sum(pr.x)
end

dg(x)=ForwardDiff.derivative(g,x)


julia> f(1.0)
2.7401598813234886

julia> df(1.0)
-2.1873699271126963

julia> g(1.0)
2.7401598813234886

julia> dg(1.0)
-2.1873699271126963
6 Likes