# LinearSolve.jl for many values of b?

**URL:** <https://discourse.julialang.org/t/linearsolve-jl-for-many-values-of-b/94980>\
**Category:** Performance\
**Tags:** linearsolve\
**Created:** [February 21, 2023, 9:40pm UTC](https://discourse.julialang.org/t/linearsolve-jl-for-many-values-of-b/94980 "2023-02-21T21:40:15Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![moble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moble/32/23535_2.png) [@moble](https://discourse.julialang.org/u/moble)\
**Post date:** [February 21, 2023, 9:40pm UTC](https://discourse.julialang.org/t/linearsolve-jl-for-many-values-of-b/94980/1 "2023-02-21T21:40:16Z")

</div>

I’m writing a package and I really like the idea of LinearSolve.jl, because I want my code to perform well and be flexible: work on CPU or GPU, be differentiable, not allocate, etc. But one of my most important use cases involves solving `A\b` for many `b` vectors. Without LinearSolve, I would just pre-compute `lu(A)`, and then farm out `ldiv!`s to many threads. I’d like to do the same with LinearSolve, but it looks like the interface forbids it — `set_b` obviously isn’t thread safe in that sense.

To be a little more concrete for my application, `A` will typically have a size in the range 70-200 (squared), and there will be ~50,000 solves each time. I use machines with 100+ cores per node, so I can get silly with the threads. But I also have to repeat the whole thing (serially) many thousands of times with different `A`s.

Is this just a case where LinearSolve was designed with very different goals in mind, making the built-in way makes more sense for my application? Maybe I should just use `RecursiveFactorization.lu!` for some Julianity? Maybe even TriangularSolve.jl?

* * *

Here’s a MWE reflecting what I had hoped to be able to do:

```julia
n=4
N=10
A = rand(n, n)
b = rand(n, N)
prob = LinearProblem(A, b)
linsolve = init(prob)
sol = solve(linsolve) # Threading somehow???

```

The last line fails with a BoundsError as `ldiv!` tries to set `n*N` elements of a vector with `n` elements. On the other hand `A\b` does what I want.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [December 16, 2024, 12:52pm UTC](https://discourse.julialang.org/t/linearsolve-jl-for-many-values-of-b/94980/2 "2024-12-16T12:52:40Z")

</div>

I’m looking for exactly the same thing. Any updates for 2024?

---

<div class="post-metadata">

**Author:** ![moble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moble/32/23535_2.png) [@moble](https://discourse.julialang.org/u/moble)\
**Post date:** [December 16, 2024, 2:20pm UTC](https://discourse.julialang.org/t/linearsolve-jl-for-many-values-of-b/94980/3 "2024-12-16T14:20:47Z")

</div>

[Looks like no](https://github.com/SciML/LinearSolve.jl/issues/552#issuecomment-2440737551), as of Oct. 28 at least — though Chris Rackauckas is aware of the need, and has a design in mind, so there’s hope.

I ended up just using the manual approach, but would still be interested in something better.