Using LinearMaps.jl Inversemap for Preconditioning IterativeSolvers.jl

General Usage
,
1

I’ve seen some examples of creating a matrix-free method via linear operators, structs, and linear maps for IterativeSolvers.jl and LinearMaps.jl . However, I haven’t seen any examples of creating matrix-free preconditioner for solving systems.

So here’s what I’m trying to do. Say I want to solve a system:

\underline{A}\cdot \bar{x} = \bar{b}

With a preconditioning matrix ‘H’ such that:

\underline{H}\cdot\bar{z} = \bar{v}

Then the iterative solver would do something like:

\bar{x}^{k+1} = \bar{x}^{k} + w^{k} \cdot \bar{z}^{k} = \bar{x}^{k} + w^{k} \cdot \underline{H}^{-1} \bar{v}^{k}

Where ‘v^{k}’ could be the residual:

\bar{v} = \bar{r}

So according to the IterativeSolvers.jl, the preconditioner needs to have a method that allos for ldiv! and P\x

So if I have a preconditioning function for \underline{H} such as:

Hpre = ( x::Array{Float64}) -> LinearMap(length(x); ismutating=true) do V,Z

#.... do matrix-free H
return Z
end

Then would it be possible to just do:

Hpresolve = (x,H,b) -> IterativeSolvers.sor!(fill!(x, 0), H, b)


Hinvpre = InverseMap(Hpre; solver=Hpresolve)

S = (x) -> LinearMap(length(x); ismutating=true) do X,B
#some linear operator S*X = B
return B
end

#rhs
b = zeros(#something)

x = IterativeSolvers.cg(S, b, Pl = Hinvpre)

I’d like to get input on this since I’m having trouble figuring out how I would implement this.

2

@erny123 I don’t know well LinearMaps.jl but you can find matrix-free preconditioners that should work with IterativeSolvers.jl here.
Note that I only tested them with Krylov.jl.

1 Like