Skip to content

verbose setting changes behavior of LU/QR switching #739

@hexaeder

Description

@hexaeder

Describe the bug 🐞

When default LU factorization fails, NonlinearSolve falls back to QR. However the behavior depends on the log level.

small systems (non blas)

  • verbose=false -> silent switching
  • verbose=true -> silent switching (used to be a warning per default)
  • verbose=Detailed() -> warning about withcin

larger sytsems (blas)

  • verbose=false -> silent switching
  • verbose=true -> fails with blas LU singularity error

Expected behavior

Don't change auto-switching behavior when changing log-level.

Minimal Reproducible Example 👇
small system singular jacobian

function nonlin(du, u, p)
    x, y = u
    du[1] = x + y - 1          # equation 1
    du[2] = 2*(x + y - 1)      # redundant equation → causes singular Jacobian
end
u0 = [0.2, 0.55]
prob = NonlinearProblem(nonlin, u0)
solve(prob, verbose=false) # works
solve(prob, verbose=true) # works, used to warn
solve(prob, verbose=Detailed()) # warns

blas sized system with singular jacobian

function nonlin_blas(du, u, p)
    f1 = sum(u) - 1
    du[1] = f1
    for i in 2:length(u)
        du[i] = i * f1       # makes Jacobian rank-1 → LU fails at solve step
    end
end
u0_blas = ones(17) .* 0.2    # initial state (nonsingular? doesn't matter: J is singular everywhere)
prob_blas = NonlinearProblem(nonlin_blas, u0_blas)
solve(prob_blas, verbose=false) # works
solve(prob_blas, verbose=true) # blas error
solve(prob_blas, verbose=Detailed()) # blas error

Error & Stacktrace ⚠️

┌ Error: BLAS/LAPACK dgetrf: Matrix is singular
│   Details: U(3,3) is exactly zero. The factorization has been completed, but U is singular and division by U will produce infinity.
│   Return code (info): 3
│   matrix_size: (17, 17)
│   rhs_type: Vector{Float64}
│   rhs_size: (17,)
│   memory_usage_MB: 0.0
│   matrix_type: Matrix{Float64}
│   element_type: Float64
└ @ LinearSolve ~/.julia/packages/LinearSolve/ZgNCX/src/mkl.jl:267
ERROR: BLAS/LAPACK dgetrf: Matrix is singular
  Details: U(3,3) is exactly zero. The factorization has been completed, but U is singular and division by U will produce infinity.
  Return code (info): 3
  matrix_size: (17, 17)
  rhs_type: Vector{Float64}
  rhs_size: (17,)
  memory_usage_MB: 0.0
  matrix_type: Matrix{Float64}
  element_type: Float64
Stacktrace:
  [1] emit_message(message::String, level::Base.CoreLogging.LogLevel, file::String, line::Int64, _module::Module)
    @ SciMLLogging ~/.julia/packages/SciMLLogging/fDS7T/src/utils.jl:57
  [2] macro expansion
    @ ~/.julia/packages/SciMLLogging/fDS7T/src/utils.jl:123 [inlined]
  [3] solve!(cache::LinearSolve.LinearCache{…}, alg::MKLLUFactorization; kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/ZgNCX/src/mkl.jl:267
  [4] solve!
    @ ~/.julia/packages/LinearSolve/ZgNCX/src/mkl.jl:239 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/LinearSolve/ZgNCX/src/default.jl:514 [inlined]
  [6] solve!(::LinearSolve.LinearCache{…}, ::LinearSolve.DefaultLinearSolver; assump::OperatorAssumptions{…}, kwargs::@Kwargs{})
    @ LinearSolve ~/.julia/packages/LinearSolve/ZgNCX/src/default.jl:503
  [7] solve!
    @ ~/.julia/packages/LinearSolve/ZgNCX/src/default.jl:503 [inlined]
  [8] #solve!#21
    @ ~/.julia/packages/LinearSolve/ZgNCX/src/common.jl:441 [inlined]
  [9] solve!
    @ ~/.julia/packages/LinearSolve/ZgNCX/src/common.jl:440 [inlined]
 [10] #_#1
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/ext/NonlinearSolveBaseLinearSolveExt.jl:24 [inlined]
 [11] LinearSolveJLCache
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/ext/NonlinearSolveBaseLinearSolveExt.jl:14 [inlined]
 [12] solve!(cache::NonlinearSolveBase.NewtonDescentCache{…}, J::Matrix{…}, fu::Vector{…}, u::Vector{…}, idx::Val{…}; skip_solve::Bool, new_jacobian::Bool, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/4j4KD/src/descent/newton.jl:119
 [13] solve! (repeats 2 times)
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/descent/newton.jl:94 [inlined]
 [14] step!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing)
    @ NonlinearSolveFirstOrder ~/.julia/packages/NonlinearSolveFirstOrder/inFag/src/solve.jl:266
 [15] step!
    @ ~/.julia/packages/NonlinearSolveFirstOrder/inFag/src/solve.jl:244 [inlined]
 [16] #step!#180
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:550 [inlined]
 [17] step!
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:544 [inlined]
 [18] solve!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:278
 [19] __solve(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:264
 [20] __solve
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:261 [inlined]
 [21] macro expansion
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:469 [inlined]
 [22] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias::SciMLBase.NonlinearAliasSpecifier, verbose::NonlinearVerbosity{…}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:414
 [23] __generated_polysolve
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:414 [inlined]
 [24] #__solve#166
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:387 [inlined]
 [25] __solve
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:384 [inlined]
 [26] #__solve#16
    @ ~/.julia/packages/NonlinearSolve/JTN4J/src/default.jl:27 [inlined]
 [27] __solve
    @ ~/.julia/packages/NonlinearSolve/JTN4J/src/default.jl:24 [inlined]
 [28] #__solve#167
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:398 [inlined]
 [29] __solve
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:390 [inlined]
 [30] #solve_call#150
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:162 [inlined]
 [31] solve_call
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:128 [inlined]
 [32] #solve_up#149
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:116 [inlined]
 [33] solve_up
    @ ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:109 [inlined]
 [34] solve(::NonlinearProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, verbose::Detailed, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/4j4KD/src/solve.jl:87
 [35] top-level scope
    @ REPL[196]:1
 [36] top-level scope
    @ REPL:1
Some type information was truncated. Use `show(err)` to see complete types.

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions