# Best performance for initialising a variable number of matrices/vectors within a Crank-Nicolson scheme

**URL:** <https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083>\
**Category:** Performance\
**Tags:** question, solver\
**Created:** [April 17, 2024, 3:47pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083 "2024-04-17T15:47:37Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![lhcgeneva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lhcgeneva/32/25595_2.png) [@lhcgeneva](https://discourse.julialang.org/u/lhcgeneva)\
**Post date:** [April 17, 2024, 3:47pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/1 "2024-04-17T15:47:37Z")

</div>

I’m solving a diffusion problem of several connected domains. I use a Crank-Nicolson scheme on each domain, and at every time step I connect the domains by some matching conditions at the boundary. To do the Crank-Nicolson I have to do u = A\b, at every time step for each domain, where u is the solution vector.

Now I’d like to make the code general in terms of the number of domains. Previously, for n=3 I used to have something like

```julia
A_1 = ...
A_2 = ...
A_3 = ...

```

similar for u and b, with the difference that u and b change at every time step.  
I’d then do

```julia
u_1 = A_1\b_1
u_2 = A_2\b_2
u_3 = A_3\b_3

```

When going to the arbitrary n case, I naively just used list comprehensions for creating all the matrices and vectors, however, that’s associated with a 10-20% performance hit. I was wondering about the best strategy here.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [April 17, 2024, 8:52pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/2 "2024-04-17T20:52:08Z")

</div>

Welcome to Julia Discourse!

You will find many experts here on numerical linear algebra and solutions of partial differential equations. I am not one of them, but I do have a couple of suggestions.

I think you might get more responses if you showed in more explicit detail what you’ve attempted so far for the general case of `n` regions.

Also, if the `A` matrices are independent of time step, then I think you could get large execution time savings by precomputing and reusing the factorizations of these matrices. I.e., if `make_region_A(i)` is called to create the `A` matrix for region `i`, then I would replace, e.g.,

```julia
As = [make_region_A(i) for i in 1:n]

```

with

```julia
factored_As = [factorize(make_region_A(i)) for i in 1:n]

```

and then for each time step the solves would look like

```julia
for i in 1:n
    ldiv!(factored_As[i], b[i])
end

```

This avoids unnecessary allocations of new `u` vectors at each time step.

I’m sure the experts will chime in with lots of additional (and better) suggestions.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 17, 2024, 8:59pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/3 "2024-04-17T20:59:56Z")

</div>

> [@lhcgeneva](#):
>
> with the difference that u and b change at every time step.  
> I’d then do
> 
> ```julia
> u_1 = A_1\b_1
> u_2 = A_2\b_2
> u_3 = A_3\b_3
> 
> ```

If the matrices stay the same, while only u and b change on each timestep, then you should definitely precompute the matrix factorizations.

e.g. precompute `F_1 = lu(A_1)` (or `cholesky` if it is PSD), and then do `u_1 = F_1 \ b_1`. Better yet, update `u` in-place with `ldiv!(u_1, F_1, b_1)`.

(Oh, I see that @PeterSimon already made a similar suggestion.)

> [@lhcgeneva](#):
>
> I naively just used list comprehensions for creating all the matrices and vectors, however, that’s associated with a 10-20% performance hit.

If you’re spending lots of time allocating vectors, then you should probably be allocating them once and then working in-place.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 17, 2024, 9:43pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/4 "2024-04-17T21:43:14Z")

</div>

Also, is there a specific reason you are using hand-written Crank-Nicolson? You likely will get better results with `OrdinaryDiffEq.jl` which would let you use much fancier solvers easily.

---

<div class="post-metadata">

**Author:** ![lhcgeneva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lhcgeneva/32/25595_2.png) [@lhcgeneva](https://discourse.julialang.org/u/lhcgeneva)\
**Post date:** [April 19, 2024, 6:56am UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/5 "2024-04-19T06:56:05Z")

</div>

Great, thanks everyone! The factorisation / lu decomposition made a 20% difference already. @Oscar_Smith : I’m interested in this! So here’s the problem stated a bit more formally. Could you give me a few keywords/hints on how to use OrdinaryDiffEq.jl effectively to solve this quicker than using backslash?

 ![Screenshot 2024-04-19 at 08.54.43](https://global.discourse-cdn.com/julialang/original/3X/5/5/554f643edbe294180901add738baad8eeee01447.png)  
 ![Screenshot 2024-04-19 at 09.01.23](https://global.discourse-cdn.com/julialang/original/3X/f/7/f75e9c95dfaca390abdd5d36dbbd31d8b1c79ccf.png)

I thought the biggest challenge was to use an efficient solver while at the same time maintaining enough flexibility to update the matching conditions between domains at every time step. Otherwise, of course, a variable time step, for example, could hugely improve things…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 19, 2024, 12:35pm UTC](https://discourse.julialang.org/t/best-performance-for-initialising-a-variable-number-of-matrices-vectors-within-a-crank-nicolson-scheme/113083/6 "2024-04-19T12:35:57Z")

</div>

5 posts were split to a new topic: [Why re-use factorization(A) rather than inv(A)](https://discourse.julialang.org/t/why-re-use-factorization-a-rather-than-inv-a/113234)