Sorry that this reply is a little late.... In addition to the other
replies, this seems like the sort of problem that might be handled by
macros: the macro can rewrite a more natural indexing strategy into the
array's index and can automatically add boilerplate that checks if the
array is the right size, etc. Writing a `@fortran_indexing` macro and
putting it in front of your first loop might be enough to convert some code
directly.
As an example, I wrote a toy macro for a presentation that does something
that might be similar: automatically generate loops to generate random
variables from an ARMA time series model. I.e., when you write
@loop_ts y[t+1] = 0.8y[t] + 0.02y[t-2] + e[t+1]
it automatically replaces it with the code
for _t in 3:(length(y) - 1)
y[_t+1] = 0.8y[_t] + 0.02y[_t-2] + e[_t+1]
end
Like I said, it's just a toy, but I've put the code in a gist if the idea
useful:
https://gist.github.com/grayclhn/d948a9d11eb05339ebaa
(I wrote this code a year ago, but it still seems to work fine.) It would
be pretty straightforward to extend the code to create they `y` vector
automatically, generate e[t+1] from a specified distribution, etc.
Macros are described in the metaprogramming section of the user's guide:
http://julia.readthedocs.org/en/latest/manual/metaprogramming/
--Gray
On Tuesday, September 29, 2015 at 2:04:30 PM UTC-5, Mark Sherlock wrote:
>
> Thanks for all the replies, I'll check the links provided.
> Just to give a summary of the reasons.... As John Gibson pointed out
> Fourier expansions would be one example.
>
> When solving PDE's we use "ghost cells" to provide boundary conditions. A
> physical variable, let's say density, could be represented in space by a 3D
> array,
> with each spatial index (ix,iy,iz) taking on values like ix=1..nx,
> iy=1..ny, iz=1..nz etc. If the boundary is, say, reflective, then you need
> to set the value of the density in
> the part of space you can't "see", e.g. in the regions ix=-2..0 and
> ix=nx+1..nx+3
>
> Obviously this can be handled by offsetting each array index by some
> amount. The problem is mainly that what we have written down on paper, in
> books and papers on
> numerical analysis etc are numerical equations which naturally count from
> 1 to n along some dimension. Couple that to the fact that many codes rely
> on arrays that
> consist of mixtures of real space (3D say) and special function
> expansions, with Fourier expansion being one example.
> I for example use a 7-dimensional array to describe some problems (3
> dimensions for real space, and 3 dimensions for velocity space). Real space
> is not expanded,
> so it counts from 1..n naturally, while velocity space has 2 dimensions
> expanded with one type of series (counting from 0..n) and the other
> dimension expanded in
> a different series (counting from -n..n). The final seventh dimension
> naturally has only 2 indices 0 and 1 (which accounts nicely for the real
> and imaginary parts of the function).
>
> Within each array loop, it is nice to be able to use if statements that
> correspond to what you have written down on paper, e.g.:
>
> if(in>0) then blah
> if(im<1) then something else
>
> or something like:
>
> do i=0,1
> do im=-mmax,mmax
> do in=max(im,0),nmax
> do ir=1,nr
> do ix=1,nx
> do iy=1,ny
> do iz=2,nz-1
>
> blah
>
> enddo
> enddo
> enddo
> enddo
> enddo
> enddo
> enddo
>
> (incidentally, imagine how many tabs would be required to write the above
> in python!).
>
> Within a code, there may be many such loops and structures with different
> looping conditions.
>
> When you have a code with 20000-40000 lines that you want to update to a
> new language, having to change all of these to account for offsetting is
> too much work for anyone to seriously consider.
> Not to mention the fact you can no longer easily code up the difference
> equations you have written on paper (and possible spent months creating).
> That in a nutshell is why nobody is switching from
> Fortran (even the young guys are writing new codes in Fortran). There are
> other nice things about Fortran arrays, such as the ability to quickly set
> all elements of an array to zero with:
>
> a=0.0
>
> I hope that sort of makes it clear the type of problems we deal with?!
>
> Cheers
>
