**Or is there something wrong with my code?**

Below is a MVE of my code. `func1`

and `func2`

do the same thing. Why is `func2`

faster than `func1`

?

Why can’t the compiler figure it out? I think `func1`

more clearly expresses my intent of performing a vectorized operation. Is there something in `func1`

that prevents the compiler from generating good code?

I am using `Julia Version 1.3.1`

, `Commit 2d57411 (2019-12-30 21:36 UTC)`

.

```
function func1(a, b, c)
for jdx = 2:size(a, 2)
a[:, jdx] .= a[:, jdx - 1] .+ b[:, jdx] .* c[jdx]
end
return a
end
function func2(a, b, c)
for jdx = 2:size(a, 2)
for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
```

With

```
a = zeros(3, 100)
b = randn(3, 100)
c = randn(100)
```

The results are:

```
@benchmark func1(a, b, c)
BenchmarkTools.Trial:
allocs estimate: 198
--------------
median time: 6.413 μs (0.00% GC)
--------------
```

and

```
@benchmark func2(a, b, c)
BenchmarkTools.Trial:
allocs estimate: 0
--------------
median time: 407.980 ns (0.00% GC)
--------------
```

For a comprehensive comparison—and an example of what a newbie tries—, see below: (PS: interpolation as in `func1($a, $b, $c)`

does not seem to matter, contrary to what we are told).

**Time in seconds**

Function | No interp. | Interpolation | Interp. vs. no-interp |
---|---|---|---|

func1 | 6.41E-06 | 6.20E-06 | 3.50% |

func2 | 4.08E-07 | 3.97E-07 | 2.79% |

func1b | 6.42E-06 | 6.19E-06 | 3.60% |

func2b | 4.08E-07 | 3.97E-07 | 2.85% |

func1c | 6.28E-06 | 6.15E-06 | 2.06% |

func2c | 3.99E-07 | 3.82E-07 | 4.53% |

func2d | 3.83E-07 |
3.72E-07 |
2.75% |

func1e | 4.84E-06 | 4.80E-06 | 0.77% |

func2e | 4.01E-07 | 3.89E-07 | 3.04% |

```
function func1(a, b, c)
for jdx = 2:size(a, 2)
a[:, jdx] .= a[:, jdx - 1] .+ b[:, jdx] .* c[jdx]
end
return a
end
function func2(a, b, c)
for jdx = 2:size(a, 2)
for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
function func1b(a::Matrix, b::Matrix, c::Vector)
for jdx = 2:size(a, 2)
a[:, jdx] .= a[:, jdx - 1] .+ b[:, jdx] .* c[jdx]
end
return a
end
function func2b(a::Matrix, b::Matrix, c::Vector)
for jdx = 2:size(a, 2)
for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
function func1c(a, b, c)
@inbounds for jdx = 2:size(a, 2)
a[:, jdx] .= a[:, jdx - 1] .+ b[:, jdx] .* c[jdx]
end
return a
end
function func2c(a, b, c)
@inbounds for jdx = 2:size(a, 2)
for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
function func2d(a, b, c)
@inbounds for jdx = 2:size(a, 2)
@simd for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
function func1e(a, b, c)
@inbounds @views for jdx = 2:size(a, 2)
a[:, jdx] .= a[:, jdx - 1] .+ b[:, jdx] .* c[jdx]
end
return a
end
function func2e(a, b, c)
@inbounds for jdx = 2:size(a, 2)
for idx in axes(a, 1)
a[idx, jdx] = a[idx, jdx - 1] + b[idx, jdx] * c[jdx]
end
end
return a
end
```