SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
244 stars 52 forks source link

Better error message when eltype(b)<:Interger #206

Open KronosTheLate opened 1 year ago

KronosTheLate commented 1 year ago

The following code errors:

julia> using LinearSolve

julia> let 
           A =  [1.0           0.0       -1.0         0.0
           0.0       -3152.28       0.0      3152.28
           0.388658      0.921382  -1.76854    -1.45868
           0.921382     -0.388658   1.45868     1.76854]
           b = [0, 0, 1, 0]
           prob = LinearProblem(A, b)
           solve(prob)
       end
ERROR: InexactError: Int64(3.8064439666352263)
Stacktrace:
  [1] Int64
    @ ./float.jl:788 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:966 [inlined]
  [4] ldiv!(A::LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}}, b::Vector{Int64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.8.2+0.x64/share/julia/stdlib/v1.8/LinearAlgebra/src/triangular.jl:1200
  [5] ldiv!(A::LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, B::Vector{Int64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.8.2+0.x64/share/julia/stdlib/v1.8/LinearAlgebra/src/lu.jl:409
  [6] _ldiv!
    @ ~/.julia/packages/LinearSolve/E7qHt/src/factorization.jl:5 [inlined]
  [7] #solve#6
    @ ~/.julia/packages/LinearSolve/E7qHt/src/factorization.jl:21 [inlined]
  [8] solve(cache::LinearSolve.LinearCache{Matrix{Float64}, Vector{Int64}, Vector{Int64}, SciMLBase.NullParameters, Nothing, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, IterativeSolvers.Identity, IterativeSolvers.Identity, Float64, Nothing}, alg::GenericLUFactorization{LinearAlgebra.RowMaximum})
    @ LinearSolve ~/.julia/packages/LinearSolve/E7qHt/src/factorization.jl:16
  [9] solve(::LinearSolve.LinearCache{Matrix{Float64}, Vector{Int64}, Vector{Int64}, SciMLBase.NullParameters, Nothing, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, IterativeSolvers.Identity, IterativeSolvers.Identity, Float64, Nothing}, ::Nothing; assumptions::LinearSolve.OperatorAssumptions{Nothing}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ LinearSolve ~/.julia/packages/LinearSolve/E7qHt/src/default.jl:145
 [10] solve
    @ ~/.julia/packages/LinearSolve/E7qHt/src/default.jl:141 [inlined]
 [11] #solve#5
    @ ~/.julia/packages/LinearSolve/E7qHt/src/common.jl:158 [inlined]
 [12] solve(::LinearSolve.LinearCache{Matrix{Float64}, Vector{Int64}, Vector{Int64}, SciMLBase.NullParameters, Nothing, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, IterativeSolvers.Identity, IterativeSolvers.Identity, Float64, Nothing})
    @ LinearSolve ~/.julia/packages/LinearSolve/E7qHt/src/common.jl:157
 [13] solve(::LinearProblem{Nothing, true, Matrix{Float64}, Vector{Int64}, SciMLBase.NullParameters, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ LinearSolve ~/.julia/packages/LinearSolve/E7qHt/src/common.jl:148
 [14] solve(::LinearProblem{Nothing, true, Matrix{Float64}, Vector{Int64}, SciMLBase.NullParameters, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}})
    @ LinearSolve ~/.julia/packages/LinearSolve/E7qHt/src/common.jl:147
 [15] top-level scope
    @ REPL[2]:8

The error can be quite intimidating, however the problem is simple and bound to arise several times. Could a better error message be provided in this case?

fredrikekre commented 1 year ago

Perhaps https://github.com/SciML/LinearSolve.jl/blob/8550cf9d09bab8ab5086e7b0e09b2dee2dd6a949/src/common.jl#L107 can be updated to

u0 = similar(b, promote_type(eltype(A), eltype(b)), size(A, 2))
KronosTheLate commented 1 year ago

Should it not always be a Float64, as the process of producing the new elements involves division? Even in the case of eltype(A) = Int64?

KronosTheLate commented 1 year ago

From what I can see from the wikipedia article on invertible matrix(full disclosure, I only looked at the first) it seems like there are metods of inversion that do in fact not include division. So if such a method is used, your proposal seems perfect.

If, however, a method involving division is used, the element-type of the output will be Float64 in all cases, right? I suppose if Float32 was ever used for both matrices, it would be nice to have that as the element type of the output, which I believe could be achieved by the verbose promote_type(typeof(one(eltype(A))/one(eltype(A))), typeof(one(eltype(b))/one(eltype(b)))). Would this verbose beast not be the right type for methods involving division, and your proposal of promote_type(eltype(A), eltype(b)) would be right for methods that do not involve division? Is such a splitting of cases based on solver even possible?

ChrisRackauckas commented 1 year ago

In theory it can also be a rational number. But your promotion there should work for that as well.

KronosTheLate commented 1 year ago

The promotion with or without division? The one without should be a nobrainer to make the default AFAICT, so that could be done immideatly. That closes the issue of original example at least. I can then make a new issue for the case of eltype(A) <: Integer && eltype(b) <: Integer.

ChrisRackauckas commented 1 year ago

I think without. If you make a PR and add your test case it should be good to go.

fredrikekre commented 1 year ago

My idea with https://github.com/SciML/LinearSolve.jl/issues/206#issuecomment-1271773381 wasn't that it would cover all cases, just catch more ones. There is no good way I think to infer the correct type anyway.