Tim - many thanks for the reply. To put this into context: I have 15y+
experience with matlab, and some limited experience with other languages
(C,C++,Java,Fortran).
Here is a code snippet that brought this up. It precomputes a lot of data
that is then used in a variety of (non-standard) ways for Tight-Binding
molecular dynamics. This is a quick and dirty first-attempt implementation
to "just get it to run".
# the following arrays are generates elsewhere, d \in \{2,3\}, N is
large, alpha real
# R : d x N x N array
# E : N x N array
# VERSION 1
for a = 1:d
hHamiltonian[a, a, :, :] = slice(hHamiltonian, a, a, :, :) -
alpha * E
for b = 1:d
hHamiltonian[a,b,:,:] = slice(hHamiltonian, a, b, :, :) +
alpha^2 * E .* slice(R, a, :, :) .* slice(R, b, :, :)
end
end
instead of what I would have liked to write:
# VERSION 2
for a = 1:d
hHamiltonian[a, a, :, :] += -alpha * E
for b = 1:d
hHamiltonian[a, b, :, :] += alpha^2 * E .* R[a, :, :] .*
R[b, :, :]
end
end
Granted, since writing the above post I read up on Comprehensions (first
time I have used them, and quite like the result)
# VERSION 3
hHamiltonian = [ - alpha*E[m,n]*del[a,b] + alpha^2 * E[m,n] * R[a,m,n]
* R[b,m,n]
for a = 1:d, b=1:d, m = 1:N, n = 1:N]
I am quite happy with this last version, for the moment at least. Some
points remain:
1. what is the performance of Comprehensions compared with vectorisation or
straight for-loops?
2. The current slicing behaviour of Julia is just unexpectedly clunky.
Whether or not VERSION 2 is "good code" to write, there are many instances
where I would have written like this without a second thought. Another
example is vectorised finite element assembly which looks very similar, but
more complex.
3. VERSION 2 is still the most natural way to write for many people who do
quick and dirty numerical experiments and don't want to think too much
about good coding practises. These are the kind of people who would prefer
the code to run for 2 days rather than 2 hours, if it means they spend 1/10
of their time coding.
Any comments will be helpful. Thanks,
- Christoph
On Tuesday, 2 September 2014 11:23:34 UTC+1, Tim Holy wrote:
>
> Your example involves two tricky issues: slice behavior and the fact that,
> despite appearances, A += b is not in-place. See issues #3424, #3217, and
> precedents they link to.
>
> I'd be interested in hearing more detail about how using slice gets nasty;
> as
> you say, from this example slice doesn't look so bad. In trying to fix
> this, we
> want to make sure we're aware of all the issues.
>
> --Tim
>
> On Monday, September 01, 2014 10:02:54 PM Christoph Ortner wrote:
> > a = rand(3,3,3,3)
> > b = rand(3,3)
> > # this works:
> > a[1,1,:,:] = slice(a,1,1,:,:)+b
> > # this does not work:
> > a[1,1,:,:] += b
> >
> > This example does not look so bad, but once you use expressive variable
> > names and more dimensions it quickly gets very nasty. Because of it, I
> am
> > doing less vectorisation than I would prefer.
> >
> > I know there is a lot of discussion on slicing on the Julia issues list,
> so
> > I did not want to post another issue there.
> >
> > Is this likely to be resolved in future releases? Are there elegant
> > alternatives?
>
>
