Hello, I’m new to Julia, with prior experience of fortran and Matlab. I have a 3d magnetohydrodynamics code which does all the time stepping + spatial derivatives using 3d matrices, with the operations carried out using element by element (dot) multiplications, additions, subtractions of the coupled density matrices, 3d momentum matrices and 3d magnetic field matrices.

So for example three lines in the code would be

fxm2[iaX,iaY,iaZ] .= vx[iaXm2,iaY,iaZ].*bbyi[iaXm2,iaY,iaZ] .- vy[iaXm2,iaY,iaZ].*bbxi[iaXm2,iaY,iaZ]

fxm1[iaX,iaY,iaZ] .= vx[iaXm1,iaY,iaZ].*bbyi[iaXm1,iaY,iaZ] .- vy[iaXm1,iaY,iaZ].*bbxi[iaXm1,iaY,iaZ]

fxp1[iaX,iaY,iaZ] .= vx[iaXp1,iaY,iaZ].*bbyi[iaXp1,iaY,iaZ] .- vy[iaXp1,iaY,iaZ].*bbxi[iaXp1,iaY,iaZ]

Here all the matrices eg vx,bbyi are size (nx,ny,nz) with nx,ny,nz typically being on the order of 30-100 (though I could go higher). The indices into the matrix iaX,iaY,iaZ are the “interior” points away from the boundaries so eg iaX = nG+1:nx-nG, where nG = 2 guard channels

So I was wondering what is the fastest way of doing this in Julia. I did a few tests using gfortran -O3 compiled code versus Matlab versus Julia. For T=100 timesteps, repeated runs of the codes show the timings are roughly

```
N=31 N=91
```

Matlab. 11 sec 6 min

Julia(Float 64) 8-10 sec. 2 min

Julia(Float 32) 7-9 sec 1.5 min

F90 (real) 1.5 sec 1 min

If you’re wondering what the translations entailed, it was literally getting those lines of code to work in the relevant programming language

eg for Julia

f. xm2[iaX,iaY,iaZ] .= vx[iaXm2,iaY,iaZ].*bbyi[iaXm2,iaY,iaZ] .- vy[iaXm2,iaY,iaZ].*bbxi[iaXm2,iaY,iaZ]

for Matlab

fxm2[iaX,iaY,iaZ] = vx[iaXm2,iaY,iaZ].*bbyi[iaXm2,iaY,iaZ] .- vy[iaXm2,iaY,iaZ].*bbxi[iaXm2,iaY,iaZ]

for F90

fxm2[iaX,iaY,iaZ] = vx[iaXm2,iaY,iaZ]*bbyi[iaXm2,iaY,iaZ] - vy[iaXm2,iaY,iaZ]*bbxi[iaXm2,iaY,iaZ]

I have zero experience in parallel programming, so these codes run only on a single processor. I was wondering if there are any suggestions for speeding up the Julia code?

Thanks