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.