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
245 stars 52 forks source link

BiCGStab Fails Silently #298

Closed packquickly closed 1 year ago

packquickly commented 1 year ago

In https://github.com/SciML/LinearSolve.jl/issues/297 I mentioned that GMRES was failing silently from stagnation. Turns out BiCGStab is also failing silently from (some form of) breakdown.

Not so surprising that this fails given how I'm (mis)using the algorithm, but the failure should not be silent. Indeed, in Krylov.jl it is not silent, as Krylov.jl returns extra information which indicates the algorithm did not converge. But when called from LinearSolve.jl this information is lost and the failure becomes silent. This is also an issue in IterativeSolvers.jl where the failure is silent (see: https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/issues/338).

Reproducing code for Krylov:

import LinearSolve: KrylovJL_BICGSTAB, LinearProblem, solve 
import LinearAlgebra: cond, norm

for problem_dim in [25, 100]
    succeeded = 0
    for i = 1:100 
        condition_number = Inf
        matrix = nothing
        while condition_number > 1000
            matrix = randn(Float64, (problem_dim, problem_dim))
            condition_number = cond(matrix)
        end

        true_vec = randn(Float64, (problem_dim,))
        b = matrix * true_vec

        prob = LinearProblem(matrix, b)

        julia_soln = solve(prob, KrylovJL_BICGSTAB())

        residual_norm = norm(julia_soln - true_vec)

        succeeded += (residual_norm < 1)  # very loose tolerance!
    end
    println("problem dim: $problem_dim. succeeded: $succeeded")
end
ChrisRackauckas commented 1 year ago

Thanks for pointing this out. Are you following this up with PR?

packquickly commented 1 year ago

No, I don't plan to.