# Define mul! for DiffEqArrayOperator

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.

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.

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.

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

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:
@ 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
``````

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