# How do I write dimensionally agnostic code?

I often find myself writing very similar code for one, two, and three dimensional versions of a given operation/algorithm. Maintaining all of these versions can become tedious. Simple code generation works fairly well, but it seems as thought there must be a better way.

Is there a relatively simple way to write an operation once and have it generalize to higher or lower dimensions?

One concrete example is: suppose I need to compute the gradient of a velocity field in spectral space. In three dimensions, the Fortran loops would look something like:

do k = 1, n
do j = 1, n
do i = 1, n
phi(i,j,k) = ddx(i)*u(i,j,k) + ddx(j)*v(i,j,k) + ddx(k)*w(i,j,k)
end do
end do
end do


where the ddx array is appropriately defined. (One could also do this with matrix multiplies.) The code for a two-dimensional flow is almost exactly the same, except: the third dimension is dropped from the loops, indexes, and number of components. Is there a better way of expressing this?

Another example is: suppose I have fluid velocities defined point-wise on a three dimensional grid. To interpolate the velocity to an arbitrary location (ie, not corresponding to grid points), one can use the one-dimensional Neville algorithm successively over all three dimensions (ie, dimensional reduction). Is there an easy way to do dimensional reduction given a one-dimensional implementation of a simple algorithm?

-

## 4 Answers

You look at how deal.II (http://www.dealii.org/) does it -- there, dimension independence lies at the very heart of the library, and is modeled as a template argument to most data types. See, for example, the dimension-agnostic Laplace solver in the step-4 tutorial program:

http://www.dealii.org/developer/doxygen/deal.II/step_4.html

See also

http://dealii.sourceforge.net/index.php/Deal.II_Questions_and_Answers#Why_use_templates_for_the_space_dimension.3F

-
 I strongly agree. I haven't found a better approach than what Deal.II is doing. They use templates in a very interesting way to get around this problem. – Eldila Mar 27 '12 at 0:46 A good resource, but pretty intimidating if you don't grock C++ templates. – meawoppl Mar 27 '12 at 0:55 @Wolfgang Bangerth: Does deal.ii define iterators using templates as well? – Matthew Emmett Mar 27 '12 at 2:52 @MatthewEmmett: Yes. – Wolfgang Bangerth Mar 27 '12 at 5:47 @meawoppl: Actually, no. I regularly teach classes on deal.II, and in the beginning simply tell students that everything that says ClassWhatever<2> is in 2d, ClassWhatever<3> is in 3d, and ClassWhatever is in dim-d. I bring the lesson on templates somewhere in week 3, and while it's likely that students don't understand how it works before that, they're fully functional using it anyway. – Wolfgang Bangerth Mar 27 '12 at 5:52
show 2 more comments

This can be accomplished in any language with the following rough mental prototype:

1. Create a list of the extents of each dimension (something like shape() in MATLAB I think)
2. Create a list of your current location in each dimension.
3. Write a loop over each dimension, containing a loop over whos size changes based on the outer loop.

From there, it is a question of battling the syntax of your certain language to keep your code nd-compliant.

Having written an n-dimensional fluid-dynamics solver, I have found that it is helpful to have a language that supports unpacking a list like object as the arguments of a function. I.e. a = (1,2,3) f(a*) -> f(1,2,3). Additionally advanced iterators (such as ndenumerate in numpy) make code an order of magnitude cleaner.

-
The Python syntax for doing this looks nice and succinct. I wonder if there is a nice way to do this with Fortran... – Matthew Emmett Mar 27 '12 at 2:54
It is a bit painful to deal with dynamic memory in Fortran. Probably my major complaint with the language. – meawoppl Mar 27 '12 at 6:09

The question highlights that most "plain" programming languages (C, Fortran, at least) do not allow you to do this cleanly. An added constraint is that you want notational convenience and good performance.

Therefore, instead of writing a dimension-specific code, consider writing a code that generates a dimension-specific code. This generator is dimension-independent, even if the compute code is not. In other words, you add a layer of reasoning between your notation and the code expressing the computation. C++ templates amount to the same thing: Upside, they're built right into the language. Downside, they're somewhat cumbersome to write. This reduces the question to how to practically realize the code generator.

OpenCL lets you do code generation at run time fairly cleanly. It also makes for a very clean split between 'outer controlling program' and 'inner loops/kernels'. The outer, generating program is far less performance-constrained, and therefore might as well be written in a comfy language, like Python. That's my hope for how PyOpenCL will get used--sorry for the renewed shameless plug.

-
 Andreas! Welcome to scicomp! Glad to have you on the site, I think you know how to get in touch with me if you have any questions. – Aron Ahmadia Mar 31 '12 at 6:09 +10000 for automatic code generation as a solution to this problem instead of C++ magic. – Jeff Apr 14 at 20:25

You can write an algorithm for a $n_1 \times n_2 \times n_3$ grid and then specialize to lower dimensions by setting some $n_j=1$.

-
 So to be dimension-independent, your code needs to be written for maxdim+1 dimensions, where maxdim is the maximum possible dimension the user could ever encounter. Let's say maxdim=100. How useful is the resulting code? – Jeff Apr 14 at 20:26