# Trying to parallelize using CUSOLVERRF.jl with @threads

**URL:** <https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457>\
**Category:** GPU\
**Tags:** multithreading, cuda, threads, cudajl\
**Created:** [August 7, 2025, 4:11pm UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457 "2025-08-07T16:11:04Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![mhunke](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@mhunke](https://discourse.julialang.org/u/mhunke)\
**Post date:** [August 7, 2025, 4:11pm UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/1 "2025-08-07T16:11:05Z")

</div>

Hi everyone,

I’m working with CUSOLVERRF.jl to solve multiple sparse linear systems on the GPU, using the high-level `RFLU` API: [GitHub - exanauts/CUSOLVERRF.jl: A Julia wrapper for cusolverRF](https://github.com/exanauts/CUSOLVERRF.jl)

All systems share the **same sparsity pattern** , but have **different matrix values**. I’m trying to run multiple solves in parallel using `Threads.@threads`. My goal is to have a code that is as fast as possible. At the moment I don’t want/can use CUDA itself.

**Typical code I’m using (works serially):**

```julia
rf = RFLU(dA; symbolic=:RF) # only once  
for i in 1:N 
    dA.nzVal .= new_values[i] # update matrix values     
    lu!(rf, dA) # update numeric factorization     
    copyto!(db, b[i]) # copy RHS     
    ldiv!(rf, db) # solve A*x = b 
end 

```

This works great **in serial**.

**What I want to do:**

Run this loop **in parallel** using `@threads` to speed up processing of many systems:

```julia
rf = RFLU(dA; symbolic=:RF) # only once  
@threads for i in 1:N     
    dA.nzVal .= new_values[i]     
    lu!(rf, dA)     
    copyto!(db, b[i])     
    ldiv!(rf, db) 
end

```

But this leads to wrong results because the `rf` object is **not thread-safe**. I can’t work with locks or similar methods to avoid data-races because that would lead to a serial execution.

**My question**  
Is there a safe and efficient way to reuse the same `RFLU` symbolic factorization across threads? With these points in mind:

- Can I clone `rf` so that each thread has its own copy of the numeric buffers, but reuses the symbolic structure? Copying does not work.
- I want to avoid calling `RFLU(dA)` in each thread, because that recomputes the symbolic factorization and is slow.

Any suggestions or workarounds would be really appreciated.

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [August 8, 2025, 1:13am UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/2 "2025-08-08T01:13:23Z")

</div>

> [@mhunke](#):
>
> All systems share the **same sparsity pattern** , but have **different matrix values**. I’m trying to run multiple solves in parallel using `Threads.@threads`

how do you multi-thread something that uses GPU? that doesn’t make sense

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [August 8, 2025, 7:24am UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/3 "2025-08-08T07:24:36Z")

</div>

You should go with Distributed.jl instead since I imagine you want to use multiple machine with gpus see [Multiple GPUs · CUDA.jl](https://cuda.juliagpu.org/stable/usage/multigpu/).  
If you do want to use multithreaded on 1 gpu (which I think will be far slower than just building the big block sparse array and use full power of your gpu on it) I think you will run out of mem really quickly (see [Multhreading & GPU memory management](https://discourse.julialang.org/t/multhreading-gpu-memory-management/116297)).

PS : Last little thing maybe you know but the code you showed isn’t thread safe

---

<div class="post-metadata">

**Author:** ![mhunke](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@mhunke](https://discourse.julialang.org/u/mhunke)\
**Post date:** [August 8, 2025, 7:25am UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/4 "2025-08-08T07:25:16Z")

</div>

I multi-thread CPU code where each thread gets its own stream, and then each thread uses the `lu!` or other cuda kernels. If the kernels don’t saturate the GPU individually (which they do not), there should be parallelization.

---

<div class="post-metadata">

**Author:** ![jpsamaroo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jpsamaroo/32/46804_2.png) [@jpsamaroo](https://discourse.julialang.org/u/jpsamaroo)\
**Post date:** [August 10, 2025, 7:56pm UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/5 "2025-08-10T19:56:06Z")

</div>

This is correct, using multiple streams will (usually) increase device utilization by providing sufficient work to engage all SMs.

To your original point, you’ll need to create `rf` on each task most likely, so that you have one per task. You can use TaskLocalValues.jl to do something like:

```julia
using TaskLocalValues

...

RF = TaskLocalValue{RFLU}(()->RFLU(dA; symbolic=:RF))
@threads for i in 1:N
    rf = RF[] # Allocates a new `RFLU` per task, as needed
    ...
end

```

---

<div class="post-metadata">

**Author:** ![mhunke](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@mhunke](https://discourse.julialang.org/u/mhunke)\
**Post date:** [August 11, 2025, 4:52pm UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/6 "2025-08-11T16:52:26Z")

</div>

After some digging in the source code of CUSOLVERRF.jl I found that they have a `RFBatchedLowLevel()` function—exactly what I need. So I can restructure my code and work with that. Even better than what I thought of before. But with the new function I get a `ReadOnlyMemoryError()` from the function call ``. Here is the code I used for testing:

```julia
using CUDA, SparseArrays, CUSOLVERRF, CUDA.CUSPARSE

A_cpu = spdiagm(0 => ones(Float64, 4))
println(A_cpu)
A_cpu = SparseMatrixCSC{Float64, Int32}(A_cpu)

dA = CuSparseMatrixCSR(copy(A_cpu))
@show typeof(dA)

try
    dA.nzVal .= dA.nzVal .+ 1.0
    println("dA.nzVal is writable!")
catch e
    println("Error writing to dA.nzVal: ", e)
end

try
    rf_batched = CUSOLVERRF.RFBatchedLowLevel(dA, 4)
    println("Success, type: ", typeof(rf_batched))
catch e
    println("Failure: ", e)
    rethrow()
end

```

Any ideas?

---

<div class="post-metadata">

**Author:** ![mhunke](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@mhunke](https://discourse.julialang.org/u/mhunke)\
**Post date:** [August 14, 2025, 7:26am UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/7 "2025-08-14T07:26:06Z")

</div>

I am now using CUDSS.jl Batched LU: [Batch API · CUDSS.jl](https://exanauts.github.io/CUDSS.jl/dev/batch/) and this works fine

---

<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:** [August 19, 2025, 9:09am UTC](https://discourse.julialang.org/t/trying-to-parallelize-using-cusolverrf-jl-with-threads/131457/8 "2025-08-19T09:09:29Z")

</div>

@mhunke  
I interfaced exactly what you need a few days ago.  
With the latest release of cuDSS we can now solve a batch of linear systems with the same sparsity pattern.

I interfaced the feature and added some examples in the documentation:

> **[Uniform batch · CUDSS.jl](https://exanauts.github.io/CUDSS.jl/dev/uniform_batch/)**
>
> Documentation for CUDSS.jl.

I worked on it for a project with MadNLP.jl but I’m glad if it can be useful to other people.