# Batched LU solves (or Factorizations) with Sparse Matrices

**URL:** <https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019>\
**Category:** GPU\
**Created:** [November 9, 2023, 10:02pm UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019 "2023-11-09T22:02:20Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![schev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schev/32/51307_2.png) [@schev](https://discourse.julialang.org/u/schev)\
**Post date:** [November 9, 2023, 10:02pm UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/1 "2023-11-09T22:02:20Z")

</div>

Let’s say I have a sparse set of matrices `A1, A2, ... An`, each of type `CuSparseMatrixCSC{Float32, Int32}` along with associated RHS vectors `b1, b2, ... bn`. Is there some watch to do a batch LU solve on the GPU with CUDA.jl?

If the matrices are dense, `getrf_batched!` works fine to get the LU factors. Is there some way to do batch solves/factorizations with sparse matrices @maleadt?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [November 9, 2023, 10:47pm UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/2 "2023-11-09T22:47:51Z")

</div>

how big are the a matrices? you might just want to treat them as dense

---

<div class="post-metadata">

**Author:** ![schev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schev/32/51307_2.png) [@schev](https://discourse.julialang.org/u/schev)\
**Post date:** [November 10, 2023, 12:02am UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/3 "2023-11-10T00:02:44Z")

</div>

The matrices are large network matrices (e.g., 10k x 10k), and they are extremely sparse (e.g., 99.99%). So leaving them dense is not an option, sadly!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 10, 2023, 5:15am UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/4 "2023-11-10T05:15:54Z")

</div>

LinearSolve.jl’s SimpleGMRES has optimizations for batched solving and is compatible with GPUs. That’s likely to be the best option here, though if you can supply a good preconditioner that would be helpful.

> **[Linear System Solvers · LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/#LinearSolve.SimpleGMRES)**
>
> Documentation for LinearSolve.jl.

---

<div class="post-metadata">

**Author:** ![schev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/schev/32/51307_2.png) [@schev](https://discourse.julialang.org/u/schev)\
**Post date:** [November 10, 2023, 4:57pm UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/5 "2023-11-10T16:57:23Z")

</div>

Thanks @ChrisRackauckas! Big fan of LinearSolve.jl. Unfortunately, I am solving linear systems within the context of a nonlinear optimization problem (interior point method), where the associated matrices sequentially have worse and worse condition numbers as convergence is approached (10^20 even). Thus, iterative methods (GMRES et al.) are well-known to cause slow convergence or divergence, so they are not a great choice here. If there is a GPU-batched LU/QR solver in LinearSolve.jl, I would love to hear about it. If not, maybe I can help develop it in the future.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 10, 2023, 9:02pm UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/6 "2023-11-10T21:02:39Z")

</div>

I see. [`MKLPardisoFactorize`](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/#LinearSolve.MKLPardisoFactorize) should do batching I think?

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [April 15, 2024, 5:37am UTC](https://discourse.julialang.org/t/batched-lu-solves-or-factorizations-with-sparse-matrices/106019/7 "2024-04-15T05:37:51Z")

</div>

@schev I suggest to try CUDSS.jl ([GitHub - exanauts/CUDSS.jl](https://github.com/exanauts/CUDSS.jl)).

It provides an interface to the new sparse linear solvers of NVIDIA.

They don’t provide a routine for solving batched sparse linear systems but you can create one big sparse block diagonal matrix diag(A1, A2, …, An) with the right-hand side [b1; b2; …; bn].