# Sparse LU factorization on GPU

**URL:** <https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012>\
**Category:** GPU\
**Tags:** linearalgebra, factorization\
**Created:** [October 30, 2024, 5:58pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012 "2024-10-30T17:58:19Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [October 30, 2024, 5:58pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/1 "2024-10-30T17:58:19Z")

</div>

Hi,

I am working with a `CuSparseMatrixCSR{Float64, Int32}` matrix `A` on the GPU to solve a linear system multiple times during the training of an ML model. I’m hoping to perform a **sparse** LU factorization (with optional reordering) of `A` and reuse it for efficient repeated solves, rather than factoring each time.

I noticed the availability of `CUDA.CUSOLVER.SparseCholesky` for Cholesky factorizations but couldn’t find an equivalent for LU. Is there a way to perform a sparse LU factorization in this context, or do you know of any workarounds or alternative approaches to achieving this? Any insights or library recommendations would be greatly appreciated.

Thank you!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [October 30, 2024, 6:10pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/2 "2024-10-30T18:10:43Z")

</div>

@amontoison this one’s for you

---

<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:** [October 30, 2024, 6:35pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/3 "2024-10-30T18:35:12Z")

</div>

What you need is [CUDSS.jl](https://github.com/exanauts/CUDSS.jl).  
I provide an [example](https://github.com/exanauts/CUDSS.jl?tab=readme-ov-file#example-1-sparse-unsymmetric-linear-system-with-one-right-hand-side) in the README to demonstrate how to reuse the analysis phase, which is the most expensive part as it is performed on the CPU (cuDSS relies on METIS).

We use it in interior point methods (MadNLP.jl) to solve multiple KKT systems with the same sparsity patterns, which makes it very efficient!

---

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [October 30, 2024, 6:50pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/4 "2024-10-30T18:50:44Z")

</div>

Thanks, that’s what I really needed. I have one question regarding the example provided in the readme. I don’t know in advance how many systems I will solve, and I also don’t know the right-hand-sides in advance. If I do a second call to the solver, as done below:

1. Are all computations on the GPU?
2. Can I use the same `x_gpu` as before? Should I allocate a new one as I did or not?

```julia
using CUDA, CUDA.CUSPARSE
using CUDSS
using SparseArrays, LinearAlgebra

T = Float64
n = 100
A_cpu = sprand(T, n, n, 0.05) + I
x_cpu = zeros(T, n)
b_cpu = rand(T, n)

A_gpu = CuSparseMatrixCSR(A_cpu)
x_gpu = CuVector(x_cpu)
b_gpu = CuVector(b_cpu)

solver = CudssSolver(A_gpu, "G", 'F')

cudss("analysis", solver, x_gpu, b_gpu)
cudss("factorization", solver, x_gpu, b_gpu)
cudss("solve", solver, x_gpu, b_gpu)

### SECOND CALL HERE
x_gpu_new = zeros(T, n)
b_gpu_new = CuVector(rand(T, n))
cudss("solve", solver, x_gpu_new, b_gpu_new)

```

---

<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:** [October 30, 2024, 7:06pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/5 "2024-10-30T19:06:27Z")

</div>

1. If you make a second call to the solver, as many operations as possible are performed on the GPU. CUDSS sometimes uses the CPU, but in the same way as the GPU in an asynchronous manner, so it’s more about being faster than anything else.

2. You can reuse the same storage `x_gpu` and `b_gpu` if you want

You can also use this high-level API if you prefer:

> **[Generic API · CUDSS.jl](https://exanauts.github.io/CUDSS.jl/dev/generic/#LU)**
>
> Documentation for CUDSS.jl.

---

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [October 30, 2024, 10:51pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/6 "2024-10-30T22:51:36Z")

</div>

Thanks, the high-level API is very useful. Just one last question: any way to use the same LU factorization to solve systems “from the right”, i.e. something like `b/A` instead of `A\b`?

---

<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:** [October 31, 2024, 12:42am UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/7 "2024-10-31T00:42:07Z")

</div>

What is the system that you want to solve with `b\A`?

I didn’t overloaded `\` in CUDSS.jl if the factorization is the argument on the right but I can, if I understand what it’s doing.

If `b` is a vector, do you want to to solve `A'y = b` or is it something different?

---

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [October 31, 2024, 6:53am UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/8 "2024-10-31T06:53:19Z")

</div>

Yes, that’s correct, I want to solve systems with `A'`, and so use a factorization of the form `A'=U'D'L'`. No need of overloading, if there is a more verbose way to do it is fine as well.

---

<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:** [October 31, 2024, 7:30am UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/9 "2024-10-31T07:30:57Z")

</div>

The library cuDSS doesn’t support yet the option to solve linear systems with `A'`.  
The option is already available in the API but not functional (`CUDSS_CONFIG_SOLVE_MODE` in [cuDSS Data Types — NVIDIA cuDSS documentation](https://docs.nvidia.com/cuda/cudss/types.html#cudssconfigparam-t)).

It will probably be available in future releases.  
If you have an example of application where you need to solve systems with both `A` and `A'`, I can send an email to the cuDSS team about that.

---

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [October 31, 2024, 8:39am UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/10 "2024-10-31T08:39:28Z")

</div>

I have a loss function of the form \ell(x) = \frac{1}{2}\|A^{-1}u(x)\|\_2^2, where u:\mathbb{R}^p \to \mathbb{R}^n is some function (in my case, a neural network, but it doesn’t matter much). Then if I want to take the gradient with respect to x of my loss function, it will be \nabla\ell(x) = \mathrm{J}u(x)^{\top} A^{-\top}A^{-1} u(x). This is why I need to solve systems with A^\top as well.

---

<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:** [October 31, 2024, 4:08pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/11 "2024-10-31T16:08:09Z")

</div>

Is the product AA^T too expensive to compute or dense in your case?  
Otherwise, you can use `cholesky` on this product.  
Solving a system with AA^T is equivalent to do a product with A^{-T} A^{-1} (if A is square and nonsingular), which is what you need in the gradient of your loss function.

---

<div class="post-metadata">

**Author:** ![IvanBioli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ivanbioli/32/213050_2.png) [@IvanBioli](https://discourse.julialang.org/u/IvanBioli)\
**Post date:** [November 2, 2024, 5:09pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/12 "2024-11-02T17:09:25Z")

</div>

Unfortunately, the product AA^{\top} it is both expensive to compute and dense in my case, plus the conditioning number would square.

---

<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:** [November 2, 2024, 5:26pm UTC](https://discourse.julialang.org/t/sparse-lu-factorization-on-gpu/122012/13 "2024-11-02T17:26:39Z")

</div>

The only alternative is to perform the two factorizations (A and A^T) in that case with the release 0.3.0 of cuDSS.