# Define mul! for DiffEqArrayOperator

**URL:** <https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911>\
**Category:** Performance\
**Tags:** performance, differentialequation\
**Created:** [March 15, 2022, 11:51am UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911 "2022-03-15T11:51:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 15, 2022, 11:51am UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911/1 "2022-03-15T11:51:45Z")

</div>

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:

```julia
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

```julia
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!

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 15, 2022, 1:43pm UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911/2 "2022-03-15T13:43:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 15, 2022, 6:30pm UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911/3 "2022-03-15T18:30:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 16, 2022, 8:07am UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911/4 "2022-03-16T08:07:48Z")

</div>

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?

```julia
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[:])`

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

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 16, 2022, 11:27am UTC](https://discourse.julialang.org/t/define-mul-for-diffeqarrayoperator/77911/5 "2022-03-16T11:27:17Z")

</div>

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