Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

oscardssmith
Copy link
Contributor

as discussed in #533

@ChrisRackauckas
Copy link
Member

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

@oscardssmith
Copy link
Contributor Author

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
Copy link
Member

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

@j-fu
Copy link
Contributor

j-fu commented Aug 25, 2024

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()

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants