Using SciMLOperators in LinearSolve

Specific Domains Numerics
1

I’m trying to use FunctionOperators from SciMLOperators, to solve a linear problem using matrix free method. But I’m getting wrong result. What am I doing wrong here? Thanks in advance.

Here is the MWE.

using SparseArrays
using LinearAlgebra
using LinearSolve
using SciMLOperators

# du = A * (B * u)
f(u, p, t) = p[1] * (p[2] * u)
F = FunctionOperator(f, zeros(n), zeros(n))

n = 100
u = rand(n)
u0 = zeros(n)
A, B = rand(n, n), rand(n, n)
ref = B \ (A \ u)
p = (A, B)

relative_residue_ref = norm(u - (A * (B * ref))) / norm(u) # 6.910495879970626e-14
relative_residue_function = norm(u - f(ref, p, Float64[])) / norm(u) # 6.910495879970626e-14

prob = LinearProblem(F, u, p=p)
linsolve = init(prob)

sol1 = solve(linsolve)
relative_residue_function_operator = norm(u - (A * (B * sol1.u))) / norm(u) # 384.0793068887865
1 Like