How do I write dimensionally agnostic code?

Question

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?

Answer 1

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.

Answer 2

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

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#why-use-templates-for-the-space-dimension

Answer 3

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.

Answer 4

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$.

Answer 5

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).

Answer 6

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.

Answer 7

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.