Create solver directly from SparseMatrixCSC

[j-fu]

Hi, when looking at the code, I figured out that the Problem struct is essentially a sparse matrix struct based on a version of the linked list structure plus a right hand side vector.

When the solver is created, this structure is read out and fed to the solver (_inmatrix).

This PR provides additional methods for creating a SpkSolver from a SparseMatrixCSC without the detour through Problem.

It is at an early stage - just the first test worked, and IMHO it makes sense to discuss how to develop this.

[codecov-commenter]

Codecov Report

Base: 91.01% // Head: 90.91% // Decreases project coverage by -0.10% :warning:

Coverage data is based on head (64f93ac) compared to base (92835b5). Patch coverage: 92.54% of modified lines in pull request are covered.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## main #9 +/- ## ========================================== - Coverage 91.01% 90.91% -0.11% ========================================== Files 13 14 +1 Lines 1436 1596 +160 ========================================== + Hits 1307 1451 +144 - Misses 129 145 +16 ``` | [Impacted Files](https://codecov.io/gh/PetrKryslUCSD/Sparspak.jl/pull/9?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl) | Coverage Δ | | |---|---|---| | [src/SparseCSCInterface/SparseCSCInterface.jl](https://codecov.io/gh/PetrKryslUCSD/Sparspak.jl/pull/9/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl#diff-c3JjL1NwYXJzZUNTQ0ludGVyZmFjZS9TcGFyc2VDU0NJbnRlcmZhY2Uuamw=) | `92.50% <92.50%> (ø)` | | | [src/SparseMethod/SpkSparseSolver.jl](https://codecov.io/gh/PetrKryslUCSD/Sparspak.jl/pull/9/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl#diff-c3JjL1NwYXJzZU1ldGhvZC9TcGtTcGFyc2VTb2x2ZXIuamw=) | `77.55% <100.00%> (ø)` | | | [src/SparseSpdMethod/SpkSymFct.jl](https://codecov.io/gh/PetrKryslUCSD/Sparspak.jl/pull/9/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl#diff-c3JjL1NwYXJzZVNwZE1ldGhvZC9TcGtTeW1GY3Quamw=) | `71.33% <0.00%> (-3.34%)` | :arrow_down: | | [src/SparseSpdMethod/SpkMMD.jl](https://codecov.io/gh/PetrKryslUCSD/Sparspak.jl/pull/9/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl#diff-c3JjL1NwYXJzZVNwZE1ldGhvZC9TcGtNTUQuamw=) | `96.57% <0.00%> (+0.38%)` | :arrow_up: | Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Petr+Krysl)

:umbrella: View full report at Codecov.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

[PetrKryslUCSD]

There already was a function for ingesting as sparse matrix, insparse!. I do not think you are using it?

[j-fu]

Indeed, but insparse! does the following steps: 1) get list of nonzero indices + values

I,J,V=findnz(spm)

This needs to create the I array which expands the colptr from having n+1 entries to I having nnz entries.

2) feed the data into Problem via inaij. This builds up a linked list sparse matrix structure within Problem.

At this point the data are ready for inmatrix! .

Now, the PR provides an inmatrix etc which takes the data directly from the SparseMatrixCSC, without creating a Problem. This should be visible in the timings (which I will provide when time permits).

During the development of the PR I tried to separate these additions from the core code you developed, via multiple dispatch. The only change there is now that the p component of SparseSolver is a Union{Problem,SparseMatrixCSC which should not create a performance hit due to union splitting.

As for the sparspaklu, sparspaklu! method, this is designed after the corresponding methods in SparseArrays.jl and KLU.jl.

[PetrKryslUCSD]

I like the design of the solvers. I wonder if we can gain substantially in reducing costs of ingesting the sparse matrix with the new functionality.

[j-fu]

You are right, for this type of test problems, the gain is not large:

module csc_msolver003
using Test
using LinearAlgebra
using SparseArrays
using Sparspak
using Sparspak.SpkSparseSolver: SparseSolver, solve!
using Sparspak.SpkProblem: insparse!, outsparse
using BenchmarkTools

function solve1(spm, b)
    T=eltype(spm)
    n=size(spm,1)
    begin
        p = Sparspak.SpkProblem.Problem(n, n, nnz(spm), zero(T))
        Sparspak.SpkProblem.insparse!(p, spm);
        Sparspak.SpkProblem.infullrhs!(p, b);
        s = SparseSolver(p)
    end
    solve!(s)
    p.x
end

function solve2(spm, b)
    slv = SparseSolver(spm)
    solve!(slv,b)
end

function solve3(spm, b)
    lu=sparspaklu(spm)
    lu\b
end

function _test(;T=Float64, n=20)
    spm = sprand(T, n, n, 1/n)
    spm = -spm - spm' + 40 * LinearAlgebra.I
    b=rand(T,n)

    @btime solve1($spm,$b)
    @btime solve2($spm,$b)
    @btime solve3($spm,$b)

    nothing
end
end

julia> csc_msolver003._test(n=10000)
  40.460 ms (201 allocations: 9.48 MiB)
  39.829 ms (177 allocations: 7.81 MiB)
  39.858 ms (177 allocations: 7.81 MiB)
  20.123 ms (73 allocations: 17.90 MiB)

Concerning the buildup of the sparse matrix via a linked list structure, as done in Problem I agree with you.

I made the ExtendableSparse.jl package, were the inaij is hidden behind a setindex! method (or an updateindex! method which is kind of missing for SparseMatrixCSC).

I recall being inspired by some early contact with the George/Liu book when I first implemented this in C++. I was always wondering why this is not common practice, and everyone suggests the sparse(I,J,A) approach.

I then would suggest to keep this PR dormant, I can organize my needs in a different way and would rather focus on an allocation hunt for the generic BLAS methods, see #10. The timings above show very nicely that the numbers of allocations in your code are independent of the problem size.

I still believe that the PR would make the package easier accessible from other Julia packages. But we can wait if this will be proven by corresponding requests by other users.

EDIT: link to #10

[rayegun]

Could this be merged @PetrKryslUCSD? It would making integration with LinearSolve.jl much easier.

[PetrKryslUCSD]

@j-fu : What do you think?

[j-fu]

I indeed would suggest to merge this PR and to make sparspaklu etc. part of the Sparspak.jl API. Currently I am using a copy of the code of this PR for my purposes. I would have no problem to have this lifted to LinearSolve.jl. However it uses lots of internals of Sparspak.jl. It IMHO would be rather tedious to expose these officially in an extended API. Moreover, any modifications of the internal structures would have to be labeled as breaking then.

This PR attempts to add a few methods which share these internals with the current API but keeps the details essentially hidden.

If this PR would be merged, I would remove that code from ExtendableSparse.jl and make it dependent on the corresponding updated Version of Sparspak.jl.

May be @wimmerer has some suggestions for making the additions to the API even more Julian ?

As I wrote this code a couple of months ago, I would need a bit of time to get deeper into it again, so any more substantial modifications/reviews from my side will take a little bit of time.

[j-fu]

Ok @PetrKryslUCSD thanks for merging this ! Before tagging a new version I would like to make a couple more checks and small changes.

[j-fu]

Ok @PetrKryslUCSD Ithink I am ok with the way it is, my downstream tests work.