JuliaLinearAlgebra / Preconditioners.jl

A few preconditioners for iterative solvers.
https://julialinearalgebra.github.io/Preconditioners.jl/
Other
50 stars 11 forks source link

Incorrect results solving with Preconditioners #6

Closed urchgene closed 5 years ago

urchgene commented 5 years ago

Hi @mohamed82008 ,

This is a very good package you are building and will be very useful with IterativeSolvers.jl

I was conducting some testing and discovered inconsistencies with the solutions of the linear systems. I may be doing this wrong but I am asking because your insights will be extremely helpful.

using LinearAlgebra, IterativeSolvers, Preconditioners

A = rand(1000,1000); A = A*A'/1000; b = rand(1000);
p = CholeskyPreconditioner(A, 2);

A\b

ldiv!(p.L, b)

cg(A, b, Pl=p, maxiter=100)

All give different solutions.

Thanks.

mohamed82008 commented 5 years ago

Hi @urchgene !

Thanks for the issue. I found an embarrassing bug in the incomplete Cholesky. It's amazing how this slipped through the tests. Anyways, the following is probably what you need to check:


using LinearAlgebra, IterativeSolvers, Preconditioners

A = rand(1000,1000); A = A*A'/1000; b = rand(1000);
p = CholeskyPreconditioner(A, 2);

c = A \ b
d = p.L' \ (p.L \ b)
e = cg(A, b, Pl=p, maxiter=100)

norm(A * c - b)
norm(A * d - b)
norm(A * e - b)
mohamed82008 commented 5 years ago

Fixed in the latest release.