I am trying to use CG from the IterativeSolvers.jl solver. I have “A” as a matrix-free operator but it is build using nested functions. Is linear map works for “A” inside of which there are several function calls? I do not find an example in the test folder where I can setup CG solver when A is a function handle.
Here is how to do it. I am assuming you are on Julia v0.6.
using IterativeSolvers
struct MyLinearMap{TM}
A::TM
end
# You will need to define these functions for your type
Base.A_mul_B!(c, L::MyLinearMap, b) = A_mul_B!(c, L.A, b)
Base.:*(L::MyLinearMap, b) = L.A*b
Base.eltype(L::MyLinearMap) = eltype(L.A)
Base.size(L::MyLinearMap, i::Int) = size(L.A, i)
# Let's test it
A = sprand(1000, 1000, 0.01)
A = A + A' + 30I
op = MyLinearMap(A)
norm(cg(op, op*ones(1000)) - ones(1000))
#2.979898645518114e-5
1 Like
thanks a lot