I agree that vectorizing is very nice in some circumstances, and it will be
great if Julia can get to the point of not introducing so many temporaries.
(Like Tobi, I am not convinced this is an easy problem to solve.)
But there are plenty of cases where loops are clearer even when vectorization
is both possible and efficient. For example, surely you can't tell me that
L = similar(A)
m,n = size(A)
L[2:m-1,2:n-1] = A[3:m,2:n-1] + A[1:m-2,2:n-1] + A[2:m-1,3:n] +
A[2:m-1,1:n-2] - 4A[2:m-1, 2:n-1]
is easier to read and understand than
L = similar(A)
for j=2:size(A,2)-1, i=2:size(A,1)-1
L[i,j] = A[i+1,j] + A[i-1,j] + A[i,j+1] + A[i,j-1] - 4A[i,j]
end
Not only does the loops-version have fewer keystrokes, it's a heck of a lot
clearer that what you're doing is computing the 2d discrete laplacian.
--Tim
On Wednesday, May 21, 2014 05:09:57 AM Andreas Lobinger wrote:
> Hello colleague,
>
> On Wednesday, May 21, 2014 7:33:01 AM UTC+2, Andreas Noack Jensen wrote:
> > Please consider b += Diagonal(c) instead of diagm. Diagonal(c) only stores
> > the diagonal elements but works like diagm(c) for matrix arithmetic
> >
> > Vectorized code is always easier to understand
> >
> >
> > That is a strong statement. I have vectorised MATLAB code with repmat and
> > meshgrids that is completely unreadable, but would be fairly easy to
> > follow
> > if written as loops. I really enjoy that I can just write a loop in Julia
> > without slow execution.
>
> If you write code in Matlab (or similar) that's using repmat (nowadays
> bsxfun) and meshgrids (and reshape) you are operating on wrong side of
> vectorising your code. Most likely (and there might be good reasons for
> that) your data is not in the 'right' shape, so you reorder more than you
> calculate. For problems like this explicit loops are often the better way
> (so straightforward in julia) OR you take some thinking, why the data needs
> to be reordered for certain operations.
>
> The right side of vectorized code is the SIMD thinking which does
> operations (some times linalg) on Matrix/Vector forms, and that in the best
> sense leads to one-liners like
>
> M += ((A*B) + (d*C) + (E .* M0))
>
> which are easy to read/write/understand.
>
> Expressions like this seems to create in julia an incredible ammount of
> intermediates and temporary allocations, while actually the output array
> (M) is already allocated and the + and * operations just needed to be
> executed with the 'right' indexing.
>
> Some people tend to solve this by writing explicit loops and create half a
> page of code for the above one-liner.
