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

Using CHOLMOD.Factor as preconditionner #340

Open
tduretz opened this issue Jul 6, 2023 · 1 comment
Open

Using CHOLMOD.Factor as preconditionner #340

tduretz opened this issue Jul 6, 2023 · 1 comment

Comments

@tduretz
Copy link

tduretz commented Jul 6, 2023

Dear all,
I have recently tried to use Cholesky factors (of simpler problem) as a preconditionner to a Krylov solve, e.g.:

PC_fact = cholesky(PC) 
gmres!(x, A, b; Pl=PC_fact)

This fails as ldiv! is not defined for CHOLMOD.Factors:

ERROR: MethodError: no method matching ldiv!(::SparseArrays.CHOLMOD.Factor{Float64},

The issue is discussed [there] in more details (https://discourse.julialang.org/t/defining-a-preconditionner-for-iterativesolvers/86977). A temporary solution (defining ldiv! for CHOLMOD.Factors manually) is presented in that post but is there now a better approach to overcomes this?
thanks!

@tduretz
Copy link
Author

tduretz commented Jan 9, 2025

Hi, I would like to know what are the steps needed to make this possible.
How can I help in making this happening?

In practice, I want to use this approach to solving non-symmetric systems, just like in the MWE below:

using SparseArrays, LinearAlgebra
let 
    using SparseArrays, LinearAlgebra
    let 
        # Non-symmetric matrix
        M = sparse([2 1 0 0; 1.5 3 1.5 0; 0 1 3 1; 0 0 1 2])
    
        # RHS
        b = ones(4)
    
        # Direct solve (likely LU)
        x = M\b
    
        # Alternative
        M_fact = cholesky(0.5*(M+M'))
    
        # Naive defect correction / would be better to use gmres!(...)
        x = zeros(4)
        r = zeros(4)
        for iter=1:20
            r .= M*x - b
            x .-= M_fact\r 
            @show norm(r)
        end
    end
end

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

No branches or pull requests

1 participant