I don't see why you find the results surprising. This sort of vectorized
operation is the sort of thing where Matlab does well (all the inner loops
are C or Fortran code), and you wouldn't expect Julia code of the same
vectorized style to be faster. (The Julia code is a bit slower, but by a
small margin, only about 20%, which is well within the margin of variation
that you would expect from different implementations of similar algorithms.)
Where Matlab falls short, and where Julia can have a huge advantage, is on
problems that don't vectorize well, so that you have to write your own
inner loops.
However, even on a vectorizable problem like this, you can still often do
far better than generic vectorized code by writing your own inner loops
that are specialized to the problem at hand. For example, the following
Julia implementation of your computation is (a) about 7x faster than your
original Julia code on my machine and (b) much more generic (it handles
arrays of any type, and the sizes are inferred from the size of the
arguments, rather than using the global variable "n" in your original code
— in this way, it is very different than a similar implementation in C or
Fortran).
function matrixLinComb2{T1,T2}(A::AbstractArray{T1,3},y::AbstractArray{T2})
m,n,p = size(A)
length(y) == p || throw(ArgumentError("y should be length $p"))
C = zeros(promote_type(T1,T2), m,n)
for k = 1:p
yk = y[k]
for j = 1:n, i = 1:m
@inbounds C[i,j] += yk * A[i,j,k]
end
end
return C
end
Of course, it's also longer the vectorized code, but that's a typical
tradeoff of optimization. You only write code like this for
performance-critical functions, and then a few straighforward extra lines
are worth it.
