# Looking for DiffEQArrayOperator replacement

**URL:** <https://discourse.julialang.org/t/looking-for-diffeqarrayoperator-replacement/128494>\
**Category:** General Usage\
**Tags:** question, package, sciml\
**Created:** [April 28, 2025, 4:17pm UTC](https://discourse.julialang.org/t/looking-for-diffeqarrayoperator-replacement/128494 "2025-04-28T16:17:46Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![cardoza2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cardoza2/32/217504_2.png) [@cardoza2](https://discourse.julialang.org/u/cardoza2)\
**Post date:** [April 28, 2025, 4:17pm UTC](https://discourse.julialang.org/t/looking-for-diffeqarrayoperator-replacement/128494/1 "2025-04-28T16:17:46Z")

</div>

Hi,

I got asked to maintain a package that someone else in my research group developed, so I apologize if this question is off topic. The package I’m working on uses SciMLBase v1, but I’m looking to update to v2.

Previously there was a struct `DiffEQArrayOperator` (see link 1) that the guy who developed the code I’m working on used to instantiate a mass matrix and pass that to an `ODEFunction`. The upgrade to v2 appears to have eliminated this operator. I assume that the replacement should be a `MatrixOperator` (I saw this post: link 2), but I’m not really sure. I tried the `MatrixOperator`, but it caused a previously stable DE to go unstable, so I must have implemented something wrong. Any pointers in the correct direction would be helpful.

Spaces are added after https and the site name, because I haven’t been able to post with links included.  
Link 1: https:// github. com/SciML/SciMLBase.jl/blob/eb9299515d9c8f50d08dc6406f2a929eb23d8a06/src/operators/basic\_operators.jl#L91  
Link 2: https:// discourse .julialang. org/t/diffeqarrayoperator-is-missing/127500/2

Here is the snippet of code that’s breaking (it can be found here: https:// github .com/byuflowlab/GXBeam.jl/blob/d9543e9130f1fd9a5be82c8ccc0b9d58c4610097/src/interfaces/diffeq.jl#L257):

```julia
else

        # system indices
        indices = SystemIndices(assembly.start, assembly.stop, static=false, expanded=false)

        # default keyword arguments
        constants = (; assembly, indices, two_dimensional, structural_damping, force_scaling,
            xpfunc, pfunc, prescribed_conditions, distributed_loads, point_masses,
            linear_velocity, angular_velocity, gravity, t=0.0)

        # set state rate vector
        du = zeros(indices.nstates) # state rate vector must be zero
        u = rand(indices.nstates) # state vector can be anything
        t = rand() # time can be anything

        # mass matrix
        TF = eltype(system)
        nx = indices.nstates
        M = zeros(TF, nx, nx)
        update_mass_matrix! = (jacob, x, p, t) -> begin
            # zero out all mass matrix entries
            M .= 0.0
            # compute mass matrix
            mass_matrix!(jacob, x, p, (; constants..., t))
            # change sign of mass matrix
            M .*= -1
            # return result
            return M
        end
        # mass_matrix = SciMLBase.DiffEqArrayOperator(M, update_func = update_mass_matrix!)
        mass_matrix = SciMLBase.MatrixOperator(M, update_func = update_mass_matrix!) #Todo: This made things unstable.... 

        # residual function
        f = (resid, u, p, t) -> dynamic_residual!(resid, du, u, p, (; constants..., t))

        # jacobian function
        update_jacobian! = (jacob, u, p, t) -> dynamic_jacobian!(jacob, du, u, p, (; constants..., t))

        # jacobian prototype
        if sparse
            jac_prototype = spzeros(TF, nx, nx)
            update_jacobian!(jac_prototype, u, p, t)
        else
            jac_prototype = nothing
        end

    end

    if isnothing(xpfunc)
        odefunc = SciMLBase.ODEFunction{true,true}(f; mass_matrix = mass_matrix,
            jac = update_jacobian!, jac_prototype = jac_prototype)
    else
        odefunc = SciMLBase.ODEFunction{true,true}(f; mass_matrix = mass_matrix)
    end

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 29, 2025, 8:58am UTC](https://discourse.julialang.org/t/looking-for-diffeqarrayoperator-replacement/128494/2 "2025-04-29T08:58:28Z")

</div>

> [@cardoza2](#):
>
> ```julia
> M = zeros(TF, nx, nx)
> update_mass_matrix! = (jacob, x, p, t) -> begin
> # zero out all mass matrix entries
> M .= 0.0
> # compute mass matrix
> mass_matrix!(jacob, x, p, (; constants..., t))
> # change sign of mass matrix
> M .*= -1
> # return result
> return M
> end
> 
> ```

This piece of code is pretty suspect. Did you mean:

```julia
        update_mass_matrix! = (jacob, x, p, t) -> begin
            # zero out all mass matrix entries
            jacob .= 0.0
            # compute mass matrix
            mass_matrix!(jacob, x, p, (; constants..., t))
            # change sign of mass matrix
            jacob .*= -1
            # return result
            return jacob
        end

```

The previous version was relying on `M` and `jacob` having the same pointer, but that won’t necessarily always be the case.

---

<div class="post-metadata">

**Author:** ![cardoza2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cardoza2/32/217504_2.png) [@cardoza2](https://discourse.julialang.org/u/cardoza2)\
**Post date:** [June 14, 2025, 7:54pm UTC](https://discourse.julialang.org/t/looking-for-diffeqarrayoperator-replacement/128494/3 "2025-06-14T19:54:00Z")

</div>

Yeah, that seems like a reasonable change to me. Here’s what I ended up changing that worked for my package:

```julia
update_mass_matrix! = (jacob, x, p, t) -> begin
            # zero out all mass matrix entries
            jacob .= 0.0 

            # compute mass matrix
            mass_matrix!(jacob, x, p, (; constants..., t))

            # change sign of mass matrix
            jacob .*= -1
 end

        M = zeros(TF, nx, nx)
        mass_matrix = SciMLBase.MatrixOperator(M, update_func! = update_mass_matrix!) 

```

Thanks Chris