# 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?

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. Mar 31, 2012 at 6:09
• +10000 for automatic code generation as a solution to this problem instead of C++ magic. Apr 14, 2013 at 20:25

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

• 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. Mar 27, 2012 at 0:46
• A good resource, but pretty intimidating if you don't grock C++ templates. Mar 27, 2012 at 0:55
• @Wolfgang Bangerth: Does deal.ii define iterators using templates as well? Mar 27, 2012 at 2:52
• @MatthewEmmett: Yes. Mar 27, 2012 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<dim> 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. Mar 27, 2012 at 5:52

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... Mar 27, 2012 at 2:54
• It is a bit painful to deal with dynamic memory in Fortran. Probably my major complaint with the language. Mar 27, 2012 at 6:09

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? Apr 14, 2013 at 20:26

The clear answers if you want to keep Fortran speed are to use a language which has proper code generation like Julia or C++. C++ templates have already been mentioned, so I'll mention Julia's tools here. Julia's generated functions let you use its metaprogramming to build functions on demand via type information. So essentially what you can do here is do

@generated function f(x)
N = ndims(x)
quote
# build the code for the function
end
end


and then you use the N to programmatically build the code you'd want to execute given that it's N dimensional. Then Julia's Cartesian library or packages like Einsum.jl expressions can be easily built for the N dimensional function.

What's nice about Julia here is that this function is statically compiled and optimized for each new dimensional array you use, so it won't compile more than you need yet it'll get you the C/Fortran speed. In the end this is similar to using C++ templates, but it's a higher level language with a lot of tools to make it easier (easy enough that this would be a nice homework problem for an undergrad).

Another language which is good for this is a Lisp like Common Lisp. It is easy to use since like Julia it gives you the compiled AST with a lot of built in introspection tools, but unlike Julia it won't automatically compile it (in most distributions).

I am in the same (Fortran) boat. Once I have my 1D, 2D, 3D and 4D (I do projective geometry) elements I create the same operators for each type and then write my logic with high level equations that make it clear what is going on. It is not as slow as you might think to have separate loops of each operation and lots of memory copy. I let the compiler/processor do the optimizations.

For example

interface operator (.x.)
module procedure cross_product_1x2
module procedure cross_product_2x1
module procedure cross_product_2x2
module procedure cross_product_3x3
end interface

subroutine cross_product_1x2(a,b,c)
real(dp), intent(in) :: a(1), b(2)
real(dp), intent(out) :: c(2)

c = [ -a(1)*b(2), a(1)*b(1) ]
end subroutine

subroutine cross_product_2x1(a,b,c)
real(dp), intent(in) :: a(2), b(1)
real(dp), intent(out) :: c(2)

c = [ a(2)*b(1), -a(1)*b(1) ]
end subroutine

subroutine cross_product_2x2(a,b,c)
real(dp), intent(in) :: a(2), b(2)
real(dp), intent(out) :: c(1)

c = [ a(1)*b(2)-a(2)*b(1) ]
end subroutine

subroutine cross_product_3x3(a,b,c)
real(dp), intent(in) :: a(3), b(3)
real(dp), intent(out) :: c(3)

c = [a(2)*b(3)-a(3)*b(2), a(3)*b(1)-a(1)*b(3), a(1)*b(2)-a(2)*b(1)]
end subroutine


To be used in equations like

m = e .x. (r .x. g)  ! m = e×(r×g)


where e and r and g can have any dimensionality that makes mathematical sense.

A meshless numerical library C++ Medusa is claimed to be dimension indenpendent. The implementation itself strongly depends on the C++ templates.

For example in this paper they use Medusa to sole a Poisson problem in 1, 2, 3, and 4 dimensional space, with basically the same code.