Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multiple right-hand sides? #31

Open
mschytt opened this issue Mar 27, 2024 · 3 comments
Open

Multiple right-hand sides? #31

mschytt opened this issue Mar 27, 2024 · 3 comments

Comments

@mschytt
Copy link

mschytt commented Mar 27, 2024

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
Copy link
Owner

I'll have a look.

@PetrKryslUCSD
Copy link
Owner

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
Copy link
Author

mschytt commented Mar 31, 2024

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```.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants