JuliaLinearAlgebra / LinearMaps.jl

A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
Other
303 stars 42 forks source link

How to use InverseMap for preconditioner #225

Open erny123 opened 7 months ago

erny123 commented 7 months ago

Are there any examples to use InverseMap as a preconditioner?

From this links:

https://discourse.julialang.org/t/how-to-pass-a-matrix-free-preconditioner-to-iterativesolver-of-krylovkit/63315/7

https://github.com/JuliaLinearAlgebra/LinearMaps.jl/issues/81

I've done something like using the diagonal as a preconditioner:

A = LinearMap(10; issymmetric=true, ismutating=true) do C,B
    C[1] = -2B[1] + B[2]
    for i in 2:length(B)-1
        C[i] = B[i-1] - 2B[i] + B[i+1]
    end
    C[end] = B[end-1] - 2B[end]
    return C
end

Pre =  (a) -> LinearMap(10; issymmetric=true, ismutating=true) do C,B
    C[1] =0.0
    for i in 2:length(B)-1
        C[i] = a[i,i]*B[i]
    end
    C[end] = 0.0
    return C
end

aa = Matrix(-2.0*I(10))

foo(x) = InverseMap(Pre(aa); solver=IterativeSolvers.bicgstabl! )*x

struct MyPreconditioner{F}
    fun::F
end
P = MyPreconditioner(foo)

using LinearAlgebra
LinearAlgebra.ldiv!(y, P::MyPreconditioner, x) = y .= P.fun(x)
LinearAlgebra.ldiv!(P::MyPreconditioner, x) = x .= P.fun(x)

b = rand(10)

sol = bicgstabl(A,b,Pl=P, log=true, max_mv_products=400)

but this says that the solution obtaininfor naN