Hi,

I work on a toy 2D CFD solver and I wonder about the **idiomatic Julian way** to compute efficiently the product for all `j`

\in `[1:ny]`

```
Pxy[:,j]=Lx*Sxy[:j]
```

where `Lx`

is a a SymTridiagonal matrix of rank `nx`

and `Pxy`

and `Sxy`

are two 2D arrays of size `(nx,ny)`

.

I wrote a loop based function `myprod_basic`

which seems to perform significantly faster than `myprod`

. In addition `myprod`

allocates memory:

```
3.858 ms (0 allocations: 0 bytes) #myprod_basic
8.958 ms (200000 allocations: 9.16 MiB) #myprod
```

Note that this snippet may eventually be part of some material for a Julia hands-on session.

Thank you for your help.

Here is the MWE:

```
using BenchmarkTools
using LinearAlgebra
function myprod_basic(Pxy,Sxy,Lx)
nx,ny=size(Pxy)
d=Lx.dv
e=Lx.ev
for it=1:1000
# Threads.@threads for j=1:ny
for j=1:ny
@inbounds Pxy[1,j]=d[1]*Sxy[1,j]+e[1]*Sxy[2,j]
for i=2:nx-1
@inbounds Pxy[i,j]=d[i]*Sxy[i,j]+e[i-1]*Sxy[i-1,j]+e[i]*Sxy[i+1,j]
end
@inbounds Pxy[nx,j]=d[nx]*Sxy[nx,j]+e[nx-1]*Sxy[nx-1,j]
end
end
end
function myprod(Pxy,Sxy,Lx)
nx,ny=size(Pxy)
for it=1:1000
@inbounds for j=1:ny
@views mul!(Pxy[:,j],Lx,Sxy[:,j])
end
end
end
function tprod(nx,ny)
dv=rand(nx)
ev=rand(nx)
Lx=SymTridiagonal(dv,ev)
Pxy=zeros(nx,ny)
Pxyref=zeros(nx,ny)
Sxy=rand(nx,ny)
pj=Pxy[:,1]
sj=Sxy[:,1]
@btime mul!($pj,$Lx,$sj)
@btime myprod_basic($Pxyref,$Sxy,$Lx)
@btime myprod($Pxy,$Sxy,$Lx)
@assert Pxyref≈Pxy
end
nx=100
ny=100
tprod(nx,ny)
```