-
Notifications
You must be signed in to change notification settings - Fork 23
Array rewrite
We would like to switch the indexing of many of the arrays in the Fortran code to take advantage of Fortran's column-order array ordering. For instance, Fortran stores a two-dimensional array
[A_{1,1}~~A_{1,2}~~A_{1,3}]
[A_{2,1}~~A_{2,2}~~A_{2,3}]
[A_{3,1}~~A_{3,2}~~A_{3,3}]
as [A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3)]
so that the first array index moves the fastest.
Currently in two space dimensions, for example, the main solution array in classic Clawpack is dimensioned as q(1-mbc:mx+mbc, 1-mbc:my+mbc, 1:meqn)
, where the dimension indexing (i,j)
is followed by the component of the system last. Because of the way Fortran stores arrays, this means that the memory stride between q(i,j,1)
and q(i,j,2)
for example is (mx+2*mbc) * (my+2*mbc)
, leading to inefficient cache performance.
The plan is to revamp all of Clawpack so that this array will be dimensioned as q(1:meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc)
, with analogous changes to the aux
arrays. New versions of PyClaw, PetClaw, and the Python visualization routines (Visclaw) have already been written in this way. Several of the 'Classic' routines and many of the Riemann solvers have been updated as well.