MethodError in LUSolver due to Float64/Float32 mismatch when using vector_type=Vector{Float32}

Description

I am trying to generate FEM snapshots and extract the solution vectors directly in Float32. However, passing vector_type=Vector{Float32} to OrderedFESpace causes a MethodError during the time integration phase (solution_snapshots).

It appears that the underlying system matrix is factored as Float64 (resulting in an UmfpackLU{Float64, Int64} object), which clashes with the Float32 state vectors during the ldiv! operation inside the LUSolver (likely because UMFPACK strictly requires Float64 or ComplexF64).

Minimal Working Example

Here is a self-contained script that reproduces the crash:

using Gridap
using GridapROMs
using GridapROMs.ParamDataStructures

# Spatial domain and reference element
domain = (0.0, 1.0)
partition = (10,)
model = CartesianDiscreteModel(domain, partition)
reffe = ReferenceFE(lagrangian, Float32, 1) # Explicitly using Float32

# FE Space with Float32 vector type
println("Creating OrderedFESpace with Float32...")
V = OrderedFESpace(model, reffe; vector_type=Vector{Float32})

# Parametric transient setup
println("Creating TransientTrialParamFESpace...")
u(σ, t) = x -> σ[1] * x[1] * t
uₚₜ(σ, t) = parameterise(u, σ, t)
U = TransientTrialParamFESpace(V, uₚₜ)

# Weak form definition
println("Defining Weak Form & Parametric Operator...")
dΩ = Measure(Triangulation(model), 2)
m(σ, t, du, v) = ∫(v * du)dΩ
a(σ, t, u, v) = ∫(v * u)dΩ
r(σ, t, u, v) = m(σ, t, ∂t(u), v) + a(σ, t, u, v)

tdomain = 0.0:0.1:0.5
D = TransientParamSpace((1.0, 2.0), tdomain)
feop = TransientLinearParamOperator(r, (a, m), D, U, V)

# Solver and initial conditions
println("Setting up Initial Conditions and Solver...")
u₀(σ) = x -> 0.0
u₀ₚ(σ) = parameterise(u₀, σ)
uh₀ₚ(σ) = interpolate_everywhere(u₀ₚ(σ), U(σ, 0.0))

slvr = ThetaMethod(LUSolver(), 0.1, 0.5)
σₒₙ = TransientRealisation(Realisation([[1.5]]), tdomain)

# Time integration (Triggering the MethodError)
println("Extracting solution snapshots...")
snapshots, _ = solution_snapshots(slvr, feop, σₒₙ, uh₀ₚ)

Error / Stacktrace

(Note the type mismatch in the first line between UmfpackLU{Float64, Int64} and Matrix{Float32})

Creating OrderedFESpace with Float32...
Creating TransientTrialParamFESpace...
Defining Weak Form & Parametric Operator...
Setting up Initial Conditions and Solver...
Extracting solution snapshots...
ERROR: LoadError: MethodError: no method matching ldiv!(::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, ::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
The function `ldiv!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ldiv!(::IterativeSolvers.Identity, ::Any)
   @ IterativeSolvers ~/.julia/packages/IterativeSolvers/bjiMt/src/common.jl:31
  ldiv!(::LinearAlgebra.Tridiagonal, ::AbstractVecOrMat)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/tridiag.jl:1004
  ldiv!(::SparseArrays.UMFPACK.UmfpackLU{T}, ::StridedVecOrMat{T}) where T<:Union{Float64, ComplexF64}
   @ SparseArrays ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/SparseArrays/src/solvers/umfpack.jl:932
  ...

Stacktrace:
  [1] ldiv!(Y::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, A::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, B::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/factorization.jl:164
  [2] solve!(x::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, ns::Gridap.Algebra.LUNumericalSetup{SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}}, b::SubArray{Float32, 1, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true})
    @ Gridap.Algebra ~/.julia/packages/Gridap/cvsif/src/Algebra/LinearSolvers.jl:274
  [3] solve!(x::ConsecutiveParamVector{Float32, Matrix{Float32}}, ns::Gridap.Algebra.LUNumericalSetup{SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}}, A::ConsecutiveParamSparseMatrixCSC{Float32, Int64, Matrix{Float32}}, b::ConsecutiveParamVector{Float32, Matrix{Float32}})
    @ GridapROMs.ParamAlgebra ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamAlgebra/ParamSolvers.jl:41
  [4] solve!(x::ConsecutiveParamVector{Float32, Matrix{Float32}}, ls::LUSolver, A::ConsecutiveParamSparseMatrixCSC{Float32, Int64, Matrix{Float32}}, b::ConsecutiveParamVector{Float32, Matrix{Float32}})
    @ GridapROMs.ParamAlgebra ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamAlgebra/ParamSolvers.jl:25
  [5] solve!(x::ConsecutiveParamVector{Float32, Matrix{Float32}}, ls::LUSolver, op::GridapROMs.ParamODEs.ParamStageOperator, cache::GridapROMs.ParamAlgebra.SystemCache)
    @ GridapROMs.ParamAlgebra ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamAlgebra/ParamSolvers.jl:69
  [6] ode_march!(statef::Tuple{ConsecutiveParamVector{Float32, Matrix{Float32}}, ConsecutiveParamVector{Float32, Matrix{Float32}}}, solver::ThetaMethod, odeop::GridapROMs.ParamSteady.GenericParamOperator{GridapROMs.ParamODEs.LinearParamODE, GridapROMs.Utils.JointDomains}, r::GridapROMs.ParamDataStructures.TransientRealisationAt{Realisation{Vector{Vector{Float64}}}, Float64}, state::Tuple{ConsecutiveParamVector{Float32, Matrix{Float32}}, ConsecutiveParamVector{Float32, Matrix{Float32}}}, odecache::Tuple{Tuple{ConsecutiveParamVector{Float32, Matrix{Float32}}, ConsecutiveParamVector{Float32, Matrix{Float32}}}, GridapROMs.ParamAlgebra.ParamCache, GridapROMs.ParamAlgebra.SystemCache})
    @ GridapROMs.ParamODEs ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamODEs/ThetaMethod.jl:64
  [7] iterate(sol::GridapROMs.ParamODEs.ODEParamSolution{ConsecutiveParamVector{Float32, Matrix{Float32}}})
    @ GridapROMs.ParamODEs ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamODEs/ODEParamSolutions.jl:23
  [8] iterate
    @ ./iterators.jl:197 [inlined]
  [9] iterate
    @ ./iterators.jl:196 [inlined]
 [10] collect_param_solutions(sol::GridapROMs.ParamODEs.ODEParamSolution{ConsecutiveParamVector{Float32, Matrix{Float32}}})
    @ GridapROMs.ParamODEs ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamODEs/ODEParamSolutions.jl:116
 [11] macro expansion
    @ ./timing.jl:697 [inlined]
 [12] collect(sol::GridapROMs.ParamODEs.ODEParamSolution{ConsecutiveParamVector{Float32, Matrix{Float32}}})
    @ GridapROMs.ParamODEs ~/.julia/packages/GridapROMs/N03YZ/src/FEM/ParamODEs/ODEParamSolutions.jl:51
 [13] solution_snapshots(fesolver::ThetaMethod, op::GridapROMs.ParamSteady.GenericParamOperator{GridapROMs.ParamODEs.LinearParamODE, GridapROMs.Utils.JointDomains}, r::GridapROMs.ParamDataStructures.GenericTransientRealisation{Realisation{Vector{Vector{Float64}}}, Float64, Vector{Float64}}, args::Function)
    @ GridapROMs.RBTransient ~/.julia/packages/GridapROMs/N03YZ/src/RB/RBTransient/RBSolvers.jl:52

Environment

Julia version: 1.12.6
Gridap version: 0.20.8
GridapROMs version: 0.1.4

Thank you for opening the issue! The problem lies within LinearAlgebra (more specifically, the Umfpack package for factorizations). Indeed, an analogous test in Gridap will still fail, meaning the problem is not inherent to GridapROMs:

using Gridap

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)
order = 2

Ω = Triangulation(model)

degree = order
dΩ = Measure(Ω,degree)

u(x) = x[1]^2 + x[2]
f(x) = - Δ(u)(x)

reffe = ReferenceFE(lagrangian,Float32,1)
V = TestFESpace(Ω,reffe,dirichlet_tags="boundary",vector_type=Vector{Float32})
U = TrialFESpace(V,u)

a(u,v) = ∫( ∇(v)⊙∇(u) )*dΩ
l(v) = ∫( v⊙f )*dΩ 

op = AffineFEOperator(a,l,U,V)
uh = solve(op)

This will result in the following error:

julia> uh = solve(op)
ERROR: MethodError: no method matching ldiv!(::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, ::Vector{Float32})
The function `ldiv!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ldiv!(::LinearAlgebra.Tridiagonal, ::AbstractVecOrMat)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/tridiag.jl:1004
  ldiv!(::SparseArrays.UMFPACK.UmfpackLU{T}, ::StridedVecOrMat{T}) where T<:Union{Float64, ComplexF64}
   @ SparseArrays ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/SparseArrays/src/solvers/umfpack.jl:932
  ldiv!(::SparseArrays.UMFPACK.UmfpackLU{Float64}, ::StridedVecOrMat{<:Complex})
   @ SparseArrays ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/SparseArrays/src/solvers/umfpack.jl:938
  ...

Stacktrace:
 [1] ldiv!(Y::Vector{Float32}, A::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, B::Vector{Float32})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/factorization.jl:164
 [2] solve!(x::Vector{Float32}, ns::Gridap.Algebra.LUNumericalSetup{SparseArrays.UMFPACK.UmfpackLU{…}}, b::Vector{Float32})
   @ Gridap.Algebra ~/.julia/packages/Gridap/2PZdB/src/Algebra/LinearSolvers.jl:274
 [3] solve!(x::Vector{…}, ls::LUSolver, op::Gridap.Algebra.AffineOperator{…}, cache::Nothing)
   @ Gridap.Algebra ~/.julia/packages/Gridap/2PZdB/src/Algebra/LinearSolvers.jl:207
 [4] solve!(x::Vector{Float32}, nls::LUSolver, op::Gridap.Algebra.AffineOperator{SparseArrays.SparseMatrixCSC{…}, Vector{…}})
   @ Gridap.Algebra ~/.julia/packages/Gridap/2PZdB/src/Algebra/NonlinearSolvers.jl:22
 [5] solve!(u::Gridap.FESpaces.SingleFieldFEFunction{…}, solver::LinearFESolver, feop::AffineFEOperator, cache::Nothing)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/FESolvers.jl:121
 [6] solve!(uh::Gridap.FESpaces.SingleFieldFEFunction{…}, solver::LinearFESolver, op::AffineFEOperator)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/FESolvers.jl:15
 [7] solve(nls::LinearFESolver, op::AffineFEOperator)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/FESolvers.jl:53
 [8] solve(op::AffineFEOperator)
   @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/FESolvers.jl:67
 [9] top-level scope
   @ ~/git_repos/GridapROMs.jl/test/RBTransient/kolmogorov.jl:83
Some type information was truncated. Use `show(err)` to see complete types.

Until LinearAlgebra exposes new factorization methods, it's necessary to stick with Float64 element types.

Navigation Menu

Repository navigation

MethodError in LUSolver due to Float64/Float32 mismatch when using vector_type=Vector{Float32} #68

Description

Minimal Working Example

Error / Stacktrace

Environment

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

MethodError in LUSolver due to Float64/Float32 mismatch when using vector_type=Vector{Float32} #68

Description

Description

Minimal Working Example

Error / Stacktrace

Environment

Activity

nichomueller commented on Aug 4, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions