I have a sparse banded matrix A and I'd like to (direct) solve Ax=b. I have about 500 vectors b, so I'd like to solve for the corresponding 500 x's. 
I'm brand new to CUDA, so I'm a little confused as to what options I have available.

cuSOLVER has a batch direct solver cuSolverSP for sparse A\_i x\_i = b\_i using QR [here](http://devblogs.nvidia.com/parallelforall/parallel-direct-solvers-with-cusolver-batched-qr/). (I'd be fine with LU too since A is decently conditioned.) However, as far as I can tell, I can't exploit the fact that all my A\_i's are the same.

Would an alternative option be to first determine a sparse LU (QR) factorization on the CPU or GPU then perform in parallel the backsubstitution (respectively, backsub and matrix mult) on the GPU? If [cusolverSp< t >csrlsvlu()](http://docs.nvidia.com/cuda/cusolver/index.html#cusolver-lt-t-gt-csrlsvlu) is for one b\_i, is there a standard way to batch perform this operation for multiple b\_i's?

Finally, since I don't have intuition for this, should I expect a speedup on a GPU for either of these options, given the necessary overhead? x has length ~10000-100000. Thanks.
['https://Stackoverflow.com/questions/30107988', 'https://Stackoverflow.com', 'https://Stackoverflow.com/users/1672126/']
> 
> is there a standard way to batch perform this operation for multiple b\_i's?
> 
> 
> 

One option is to use the batched refactorization module in CUDA's cuSOLVER, but I am not sure if it is *standard*.

Batched refactorization module in cuSOLVER provides an efficient method to solve batches of linear systems with fixed left-hand side sparse matrix (or matrices with fixed sparsity pattern but varying coefficients) and varying right-hand sides, based on LU decomposition. Only some partially completed code snippets can be found in the official documentation (as of CUDA 10.1) that are related to it. A complete example can be found [here](https://github.com/zishun/cuSolverRf-batch).
If you don't mind going with an open-source library, you could also check out CUSP:
[CUSP Quick Start Page](https://code.google.com/p/cusp-library/wiki/QuickStartGuide)

It has a fairly decent suite of solvers, including a few preconditioned methods:
[CUSP Preconditioner Examples](http://code.google.com/p/cusp-library/source/browse/#hg%2Fexamples%2FPreconditioners)

The smoothed aggregation preconditioner (a variant of algebraic multigrid) seems to work very well as long as your GPU has enough onboard memory for it.