PetrKryslUCSD / Sparspak.jl

Direct solution of large sparse systems of linear algebraic equations in pure Julia
MIT License
36 stars 7 forks source link

Multiple right-hand sides? #31

Open mschytt opened 5 months ago

mschytt commented 5 months ago

Hi!

My goal is to use the SparseCSCInterface to get LU factorizations that solve both systems with single and multiple right-hand sides.

The documentation claims that the API supports multiple right-hand sides. The Problem interface does not complain when a matrix right-hand side is passed with infullrhs!. However, only the first column is available and used.

using SparseArrays
using Sparspak.Problem: Problem, insparse!, outsparse, infullrhs!

n = 1357
A = sprand(n, n, 1/n)
A = A + I
B = randn(n, n)

p = Problem(n, n)
insparse!(p, A)
infullrhs!(p, B)
size(p.rhs) # (1357,)

Sure enough, when solving the underlying triangularsolve! requires an AbstractVector.

Is this intended?

PetrKryslUCSD commented 5 months ago

I'll have a look.

PetrKryslUCSD commented 5 months ago

The right hand side is limited to a single vector even in the original fortran source. It would be possible to extend it to multiple righthand sides (a rectangular matrix on the right hand side). Or, one can set the righthand side vector multiple times and use a preexisting factorization for quick solves.

mschytt commented 5 months ago

I got the following methods to work for my benchmarks:

function Sparspak_LU(A_sparse)
    return Sparspak.sparspaklu(A_sparse)
end

function Sparspak_solve!(LU, b)
    # If b is a matrix loop over columns
    if size(b, 2) > 1
        for i in axes(b, 2)
            Sparspak.SparseSolver.triangularsolve!(LU, @view b[:,i])
        end
        return
    end
    return Sparspak.SparseSolver.triangularsolve!(LU, b)
end

On another note, I suppose the performance of the backsolve is up to the implementations in ``GenericBlasLapackFragments.jl```.