# Type stability of sparse solves?

**URL:** <https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235>\
**Category:** New to Julia\
**Created:** [October 17, 2025, 12:20pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235 "2025-10-17T12:20:06Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![loisel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/loisel/32/50626_2.png) [@loisel](https://discourse.julialang.org/u/loisel)\
**Post date:** [October 17, 2025, 12:20pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/1 "2025-10-17T12:20:06Z")

</div>

Hi,

If you run the following code:

```julia-auto
using SparseArrays
A = sparse([1.0f0 2.0f0; 3.0f0 4.0f0])
b = Float32[1, 2]
println(eltype(A \ b))

```

it prints `Float64`, I thought it should’ve produced a vector of `Float32`. Is this expected/known/documented?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 17, 2025, 7:57pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/2 "2025-10-17T19:57:17Z")

</div>

Note that this is not a type instability (the output type is always Float64), it is an unexpected promotion.

I never noticed this, but presumably it is because the underlying SuiteSparse library is only double precision? This issue seems relevant: [\ does not work for sparse matrices with Float32 elements · Issue #111 · JuliaSparse/SparseArrays.jl · GitHub](https://github.com/JuliaSparse/SparseArrays.jl/issues/111) … single precision Cholesky was only recently added ([Add Support for Single Precision by rayegun · Pull Request #487 · JuliaSparse/SparseArrays.jl · GitHub](https://github.com/JuliaSparse/SparseArrays.jl/pull/487)); I’m not sure about LU.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 18, 2025, 8:16pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/3 "2025-10-18T20:16:32Z")

</div>

Yup, the LU factorization is explicitly converting to `Float64` on [these lines of `SparseArrays/src/solvers/umfpack.jl`](https://github.com/JuliaSparse/SparseArrays.jl/blob/26c80c8b45dc2dca92788332a40a99b6c360d05a/src/solvers/umfpack.jl#L377-L379).

In 2023, the [SuiteSparse author wrote](https://github.com/DrTimothyAldenDavis/SuiteSparse/issues/447#issuecomment-1831536617):

> UMFPACK and SPQR are not yet revised to handle single precisions matrices

so it’s (still, apparently) a limitation of the underlying libraries.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [October 19, 2025, 2:36pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/4 "2025-10-19T14:36:39Z")

</div>

You should be able to use SparsPak w/o the implicit promotion.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 20, 2025, 2:06pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/5 "2025-10-20T14:06:01Z")

</div>

It looks like Pardiso.jl and MUMPS.jl also support single precision.

---

<div class="post-metadata">

**Author:** ![loisel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/loisel/32/50626_2.png) [@loisel](https://discourse.julialang.org/u/loisel)\
**Post date:** [November 3, 2025, 7:10pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/6 "2025-11-03T19:10:47Z")

</div>

Well, I don’t know if it’s a bug, but the way it is now just seems like a method was accidentally added for A\b with Float32 types… In my admittedly unexpert opinion, it seems like it would have been more logical not to have it at all (I think it used to be it wasn’t there), or if you’re going to silently cast it to Float64, maybe cast it back to Float32 when you return?

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [November 3, 2025, 7:42pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/7 "2025-11-03T19:42:50Z")

</div>

You can end with:

```julia-auto
julia> Float32.(A \ b) # I believe this cast rounds, not just truncates, what I want, but couldn't easily confirm, truncating might even be good enough for you?!
2-element Vector{Float32}:
 0.0
 0.5

```

While I believe you want it the operator, i.e that cast added automatically, I’m not 100% sure it’s wanted there (if the library/operator would ever returnFloat32, that will be a no-op). I.e. I suppose you get rounding to Float64, and then with the cast rounding on top of rounding (and smaller range, could it ever change to Infs or would that not happen since the input was Float32).

Double rounding is almost correct, I believe only 1, or 0, ULP off, but maybe that’s not something we want to do in a library/operator automatically. You might be perfectly ok with that though.

[I’m not expert here, I was just reading about correctly rounded math libraries, and some can do this double rounding safely, most, i.e. older do not. What I was reading didn’t have anything to do with sparse libraries, directly, but mostly likely translated to it. Julia could use such correctly rounded Float32 libraries, I think don’t and doesn’t promise it; the table makers dilemma by now only applies to larger then 32-bit, e.g. Float64, and intractable to fix, there only, my understanding.]

---

<div class="post-metadata">

**Author:** ![araujoms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/araujoms/32/217734_2.png) [@araujoms](https://discourse.julialang.org/u/araujoms)\
**Post date:** [November 5, 2025, 5:37pm UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/8 "2025-11-05T17:37:17Z")

</div>

This is an actual type instability, because if you do

```julia
using SparseArrays
A = sparse([1.0f0 2.0f0 3.0f0; 4.0f0 5.0f0 6.0f0])
b = Float32[1, 2]
println(eltype(A \ b))

```

you get `Float32`, i.e., you indeed get output of different types for inputs of the same type.

The reason is that for square `A` left division will use `lu`, which promotes to `Float64`, whereas for non-square `A` left division uses `qr`, which can handle `Float32`.

Now the promotion to `Float64` was clearly done intentionally, and it cannot be removed because that would be breaking. However, casting it back to `Float32` can be done (and should, in my opinion). I encourage you to open an issue in `SparseArrays`.

---

<div class="post-metadata">

**Author:** ![araujoms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/araujoms/32/217734_2.png) [@araujoms](https://discourse.julialang.org/u/araujoms)\
**Post date:** [April 3, 2026, 9:44am UTC](https://discourse.julialang.org/t/type-stability-of-sparse-solves/133235/9 "2026-04-03T09:44:02Z")

</div>

[Fixed now.](https://github.com/JuliaSparse/SparseArrays.jl/pull/684)