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

remove the lazy init and set isfresh-false on init #535

Open oscardssmith opened 2 weeks ago

oscardssmith commented 2 weeks ago

as discussed in https://github.com/SciML/LinearSolve.jl/pull/533

ChrisRackauckas commented 2 weeks ago

When initing you might not have real values so you might get errors due to singularities and such

oscardssmith commented 2 weeks ago

Given that the fallbacks have to factorize anyway, this doesn't seem like a problem to me. If the factorizations don't have a check=false, IMO that's a bug in the factorization.

ChrisRackauckas commented 2 weeks ago

The point of the ArrayInterface calls is to avoid the factorization and also avoid the issues with factorization failing.

j-fu commented 2 weeks ago

Hi, below is an MWE which I am interested to see working (the real application case is Newton's method with Krylov linear solvers which update preconditioners not in every step. Not sure where to put it in the moment.

using SparseArrays, LinearSolve, LinearAlgebra

struct JacobiInstance{T} 
    invdiag::Vector{T}
end

jacobi(A::AbstractMatrix{T}) where T =JacobiInstance([ 1/A[i,i] for i=1:size(A,2)])

LinearAlgebra.ldiv!(u, M::JacobiInstance, v) =  u.=M.invdiag.*v

num_prec_calls::Int=0

struct JacobiPreconditioner end

function (prec::JacobiPreconditioner)(A, p)
    global num_prec_calls
    num_prec_calls+=1
    (prec(A), I)
end

function test_jacobi(;n=10,linsolver = KrylovJL_CG(precs=JacobiPreconditioner()))
    global num_prec_calls

    num_prec_calls=0

    A=spdiagm( -1 => -ones(n-1), 0 => fill(10.0,n), 1=>-ones(n-1))
    b=rand(n)

    p = LinearProblem(A, b)
    x0=solve(p,linsolver)
    cache = x0.cache
    x0 = copy(x0)
    for i = 4:(n - 3)
        A[i, i + 3] -= 1.0e-4
        A[i - 3, i] -= 1.0e-4
    end
    LinearSolve.reinit!(cache;A)
    x1 = copy(solve!(cache, reuse_precs=true))

    # This must be true
    all(x0 .< x1) && num_prec_calls==1
end

@show test_jacobi()