# Performing mul!() with Sparse Matrices inside DiffEq.jl model

**URL:** https://discourse.julialang.org/t/performing-mul-with-sparse-matrices-inside-diffeq-jl-model/89963
**Category:** Modelling & Simulations
**Created:** [November 9, 2022, 2:24am UTC](https://discourse.julialang.org/t/performing-mul-with-sparse-matrices-inside-diffeq-jl-model/89963 "2022-11-09T02:24:54Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![jarias9](https://avatars.discourse-cdn.com/v4/letter/j/cc9497/32.png) [@jarias9](https://discourse.julialang.org/u/jarias9)
#### Post date: [November 9, 2022, 2:24am UTC](https://discourse.julialang.org/t/performing-mul-with-sparse-matrices-inside-diffeq-jl-model/89963/1 "2022-11-09T02:24:54Z")

</div>

Hello all,

I was trying to improve the runtime of a model I am running in DiffEq.jl  
I have been going through the optimization steps for an `eqnsystem!(du, u, p, t)` however I am a bit stuck with regards to how to optimize sparse matrix multiplication.

In essence I need to evaluate a flux operator at every time step and with it perform  
(Dx \* (D \* Dx)) \* u  
Here my Dx would be a sparse differentiation matrix, D is a diagonal matrix (that also updates at every timestep), and u is the solution vector.  
I want to preallocate a term Dx \* (D \* Dx) = P, where P will have a fixed sparsity pattern in the hope of reducing allocations, as it causes my code to bog down substantially.

I’ve tried things such as `mul!(P, Dx, Dx)`, but this method is extremely slow.

Currently the best performing solution is just evaluating it in-place with `P .= Dx * D * Dx` however this still creates a lot of allocations which increases substantially with the number of points. Any advice would help!

Thanks

---

<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: [November 9, 2022, 10:27am UTC](https://discourse.julialang.org/t/performing-mul-with-sparse-matrices-inside-diffeq-jl-model/89963/2 "2022-11-09T10:27:17Z")

</div>

Depending on the size, you can `modelingtoolkitize` to trace out the sparse structure down to scalarized flat equations and remove all of the sparsity overhead, and that can sometimes speed things up.

> [@jarias9](#):
>
> I want to preallocate a term Dx \* (D \* Dx) = P, where P will have a fixed sparsity pattern in the hope of reducing allocations, as it causes my code to bog down substantially.

How are you preallocating it? Is it type stable?

---

<div class="post-metadata">

### Author: ![jarias9](https://avatars.discourse-cdn.com/v4/letter/j/cc9497/32.png) [@jarias9](https://discourse.julialang.org/u/jarias9)
#### Post date: [November 9, 2022, 4:04pm UTC](https://discourse.julialang.org/t/performing-mul-with-sparse-matrices-inside-diffeq-jl-model/89963/3 "2022-11-09T16:04:50Z")

</div>

Hi @ChrisRackauckas

I had a moment of clarity and realized I didn’t have to preallocate the operator like that at all and instead could just multiply u by each operator separately (very silly of me to forget basic linear algebra). Doing this with `mul!(flux_x, Dx, u)` and just applying each operator one by one allows me to the matrix operations in place. I got around a 5-10x speed up and reduced the allocations from 100 MiB to around 1 MiB, when solving for a system with 25k unknowns.

Just to answer your questions though, I am preallocating everything by using a closure. So for example `eqnclosure!(du, u, p, t, Dx, flux_x, ...) -> eqnsystem!(du, u, p, t)` after having preallocated Dx, flux\_x, etc. Also yes everything is type stable.

Thanks for the suggestion though. I’ll keep the `modelingtoolkitize` method for tracing sparse matrices in my back pocket for a rainy day.