Define mul! for DiffEqArrayOperator

General Usage Performance
,
1

Hi,
I’m trying to use the SplitODEProblem interface in order to improve the performance of my differential equation solver.

The problem has a linear part which can sum into matrix multiplication, and a non linear part.
The variable u is multidimensional (2 dimension for now), and the matrix A is such that it needs to act on the linear index of u. Like that:

u = zeros(ComplexF64, 4,4)
A = rand(16,16)
du = A*u[:]

I want to define a DiffEqArrayOperator that will act on u in the right way which I defined.
But when the solver tries to do A*u it raises an error that the dimensions do not match.

So, I tried to define

LinearAlgebra.mul!(du::AbstractArray, L::DiffEqArrayOperator{ComplexF64, Matrix{ComplexF64}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, x::Matrix{ComplexF64}) = mul!(du, L.A, x[:])

and I still get errors. See another post below.

Any insight will be helpful.
Needless to say that changing the structure of u is not that simple.
And I tried to use LinearMaps instead of DiffEqArrayOperator with no success.
Thanks in advance!

2

I would use reshape(du, :) instead of du[:]. The latter might create a completely new copy of du so the results do not end up in the correct memory.

3

Good point!
But this would not solve my problem, would it?

4

Errors that I get when I tried different things.

  1. This is probably because du is still in the original form of u. How can I fix that?
LoadError: DimensionMismatch("A has size (16,16), B has size (16,1), C has size (4, 4)")
  1. I tried to define LinearAlgebra.mul!(...) = mul!(du[:], L.A, x[:]) 
LoadError: DimensionMismatch("A has dimensions (16,16) but B has dimensions (4,4)")
Stacktrace:
  [1] gemm_wrapper!(C::Matrix{ComplexF64}, tA::Char, tB::Char, A::Matrix{ComplexF64}, B::Matrix{ComplexF64}, _add::LinearAlgebra.MulAddMul{true, true, Bool, Bool})
    @ LinearAlgebra /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:643
  [2] mul!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:169 [inlined]
  [3] mul!
    @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/matmul.jl:275 [inlined]
  [4] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{ETDRK4{Val{:forward}, true}, true, Matrix{ComplexF64}, No......................
  1. When I tried to create MWE I got a different error
    MWE:
using DifferentialEquations, LinearAlgebra
function non_linear_part(dpsi::Array{ComplexF64,2}, psi::Array{ComplexF64,2}, t::Float64)
end

LinearAlgebra.mul!(du::AbstractArray, L::DiffEqArrayOperator{ComplexF64, Matrix{ComplexF64}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, x::Matrix{ComplexF64}) = mul!(du, L.A, x[:])

lat_size = (4,4)
psi0 = [zeros(ComplexF64,lat_size)]
tspan = (0.0, 0.1)
A = rand(16,16)
linear_part = DiffEqArrayOperator( A )

prob = SplitODEProblem(linear_part, non_linear_part, psi0, tspan);

integrator = init(prob, ETDRK4(), dt = 0.01)
sol = solve!(integrator)

Error:

LoadError: MethodError: no method matching zero(::Type{Matrix{ComplexF64}})
Closest candidates are:
  zero(::Union{Type{P}, P}) where P<:Dates.Period at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:53
  zero(::ForwardDiff.Dual) at /home/elyco/.julia/packages/ForwardDiff/CkdHU/src/dual.jl:307
  zero(::SymbolicUtils.Symbolic) at /home/elyco/.julia/packages/SymbolicUtils/lCYjx/src/types.jl:99
  ...
Stacktrace:
 [1] zero(x::Vector{Matrix{ComplexF64}})
   @ Base ./abstractarray.jl:1085
 [2] promote_f(f::SplitFunction{false, ODEFunction{true, DiffEqArrayOperator{Float64, Matrix{Float64}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, typeof(non_linear_part), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, u0::Vector{Matrix{ComplexF64}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/0PaUK/src/solve.jl:207
5

I would recommend reshaping u into a Vector from the very beginning. Then you wouldn’t need to redefine mul! at all.