Update csc interface

[PetrKryslUCSD]

@j-fu I would like to enable \ on the solver without a prior factorization. For this purpose I have replaced the logic in ldiv! with solve!. All tests pass. Please tell me what you think.

[codecov-commenter]

Codecov Report

Patch coverage: 100.00% and project coverage change: -0.06 ⚠️

Comparison is base (020b13b) 90.27% compared to head (eb80118) 90.21%.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #23      +/-   ##
==========================================
- Coverage   90.27%   90.21%   -0.06%     
==========================================
  Files          15       15              
  Lines        1831     1830       -1     
==========================================
- Hits         1653     1651       -2     
- Misses        178      179       +1     

Impacted Files	Coverage Δ
src/Problem/SpkProblem.jl	77.77% <100.00%> (+0.22%)	⬆️
src/SparseCSCInterface/SparseCSCInterface.jl	95.06% <100.00%> (+2.94%)	⬆️

... and 1 file with indirect coverage changes

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[j-fu]

But doesn't solve! redo all factorization steps ? The assumption I used is that sparspaklu returns some factorization object, and ldiv! solves the factorized problem. Due to the lack of a factorization type I (mis)used the SparseSolver itself in the corresponding state as "factorization" in the context of sparspaklu, sparspaklu! and ldiv! . LinearSolve.jl uses these things in this way, and this corresponds to what we have with e.g. with KLU.jl (which I took as an example).

I think this can be corrected by introducing a factorization type which can be returned by sparspaklu:

struct SparspakLUFactorization{IT,FT}
 s::SparseSolver{IT,FT}
end

than we are free to dispatch ldiv! either on the solver or on the factorization.

I think this is a discussion on the architecture of the package in #24 .

[PetrKryslUCSD]

ut doesn't solve! redo all factorization steps ?

It shouldn't: the sparse solver maintains a state. If something invalidates the factorization, it will factorize; otherwise not.

[j-fu]

It shouldn't: the sparse solver maintains a state. If something invalidates the factorization, it will factorize; otherwise not.

Ah ok - in that case I tinkk there would'nt be a problem. But I indeed think that we should have this factorization struct as it may be more clear in its functionality.

In the moment I am in holidays and the weather is good, so it may take some time to get to this.

[PetrKryslUCSD]

I have not discovered the reason we store the sparse matrix in the sparse solver. Do you remember what it is? We don't seem to use it once it is stored...

[j-fu]

It is for dispatching this call:

Sparspak.jl/src/SparseMethod/SpkSparseSolver.jl

Line 189 in 020b13b

success = _inmatrix!(s.slvr, s.p)

properly. If s.p is a sparse matrix, the _inmatrix! method from SparseCSCInterface is called.

[PetrKryslUCSD]

Oh, right. So, how about we embed the possibility of the discrete problem being defined by a csc matrix in the Problem? Then neither sparse LU solver, nor the LDLT SPD sparse solver need to know about the csc matrix.. Petr Krysl Prof. and Vice chair for undergraduate education Structural engineering department University of California, San Diego

…

On Wed, Mar 29, 2023 at 12:40 AM Jürgen Fuhrmann ***@***.***> wrote: It is for dispatching this call: https://github.com/PetrKryslUCSD/Sparspak.jl/blob/020b13b3a5b5e83881fda88056e581a715ea17fc/src/SparseMethod/SpkSparseSolver.jl#L189 <https://urldefense.com/v3/__https://github.com/PetrKryslUCSD/Sparspak.jl/blob/020b13b3a5b5e83881fda88056e581a715ea17fc/src/SparseMethod/SpkSparseSolver.jl*L189__;Iw!!Mih3wA!Ayk2fxL_smRzSz4qWmjmOetd_qr22CBWe9yD7vLIY3kjqNYwnm2usz2OYVDhRjK7DiZYVGecPCvlAFaksafa4Z49$> properly. If s.p is a sparse matrix, the _inmatrix! method from SparseCSCInterface is called. — Reply to this email directly, view it on GitHub <https://urldefense.com/v3/__https://github.com/PetrKryslUCSD/Sparspak.jl/pull/23*issuecomment-1488094337__;Iw!!Mih3wA!Ayk2fxL_smRzSz4qWmjmOetd_qr22CBWe9yD7vLIY3kjqNYwnm2usz2OYVDhRjK7DiZYVGecPCvlAFaksXz9kwQC$>, or unsubscribe <https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/ACLGGWBVXBEVAL7O5NNZKKDW6PRP7ANCNFSM6AAAAAAWJOBNWA__;!!Mih3wA!Ayk2fxL_smRzSz4qWmjmOetd_qr22CBWe9yD7vLIY3kjqNYwnm2usz2OYVDhRjK7DiZYVGecPCvlAFaksb2T2xHl$> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[j-fu]

So, how about we embed the possibility of the discrete problem being
defined by a csc matrix in the Problem?

I think this would be not very logical, as IMHO the Problem itself is a (linked list) sparse matrix structure together with a right hand side and a solution vector. We would have to keep these vectors uninitialized or unused when runnin with the sparselu logic.

As I see, the SparseSPDSolver in the develop branch has no p field yet but will need one in order to call the inmatrix method.

We just could do the same as for the SparseSolver and let p a Union field.
Alternatively, we could parametrize the solvers by the type of the p field.

I am ready to care about the SPD stuff for SparseMatrixCSC.

[PetrKryslUCSD]

As I see, the SparseSPDSolver in the develop branch has no p field yet but will need one in order to call the inmatrix method.

Yes. The the whole thing is half-baked at the moment.

We just could do the same as for the SparseSolver and let p a Union field.

I guess we could do that. I was thinking that perhaps code duplication could be avoided by
going through the Sparspak Problem, but I haven't thought it through entirely.

I am ready to care about the SPD stuff for SparseMatrixCSC.

That is super. Beware: the code has not been tested (very much) yet.

[j-fu]

A critical occurence of code doubling is between

Sparspak.jl/src/SparseCSCInterface/SparseCSCInterface.jl

Line 120 in 020b13b

if (inew >= xsuper[snode[jnew]])

and

Sparspak.jl/src/SparseMethod/SpkSparseBase.jl

Line 323 in 020b13b

if (inew >= xsuper[snode[jnew]])

This is the code inserting an entry into the solver data structure. This could be made into a function which can be called within both of the doit implementations. I have this change on my mental list and could make a PR for this (after the SPD stuff is stable...). However this change could incur some function call overhead so it may make sense to benchmark bevore releasing.

I think we will see a similar situation when making the CSC equivalent for

Sparspak.jl/src/SparseSpdMethod/SpkSparseSpdBase.jl

Line 233 in be5e87f

function _inmatrix!(s::_SparseSpdBase{IT, FT}, p::Problem{IT, FT}) where {IT, FT}

I don't see any other critical occurences of code doubling, though.

Once you have a running version for the SPD case with Problem I could add the CSC counterpart. I would name the factorization method sparspakldlt . I also think about changing the return value to a struct type as described above. The constructor could check the state of the solver to be factorization_done.