Which reminds me of the other generic image package that you hear about these days -- Adobe's GIL (Generic Image Library), which also support arbitrary N-dimensional storage:
http://opensource.adobe.com/gil/

Indeed the cisstVector seems nice for a user that needs all of that, but it is not what I want to do.
I working on a class/template for uniformly spaced (normally dense) multidimensional arrays with a dimension that although constant for each object once created it is not known at compile time.
The rank of the object on the other hand must be known at compile time.

Having multidimensional arrays with constant sizes (and very often you want only one or two dimensions constant) immediately increases the the complexity very much, and makes a clean design and uniform handling much less difficult (it is easy to have holes in the implementation of given corner cases).

This won't be the right class for 3D vectors, but for array of 3D vectors it might well be.

It's worth knowing that the next DMD2 release (probably in a couple of days -- I have a pre-release already) contains array operations.
So you can write stuff like:

double [] a;
double [] b;

a += b[3..$-8]*3.5;

and have it converted into a DAXPY call.
Which will mean we'll pretty much have BLAS1 built into the language.
The implications for a matrix library are quite interesting...

You mean it gets converted to something that's like a DAXPY call? It doesn't actually call a BLAS DAXPY does it? Is it configurable? Can you make it call your own routines? Or is it yet another of these "works for builtins, but not for you" things like returning lvalues from arrays.

You mean it gets converted to something that's like a DAXPY call? It doesn't actually call a BLAS DAXPY does it? Is it configurable? Can you make it call your own routines? Or is it yet another of these "works for builtins, but not for you" things like returning lvalues from arrays.

It's like a DAXPY call. It calls a new routine in the runtime (which could be forwarded to a BLAS routine if desired). Only works for built-ins. It's the first situation in which DMD will use SSE operations.

Essentially, it's just syntax sugar for BLAS1 operations: you can use an array instead of a vector class.

I think that the design is quite sound (3. iteration) and the performance reasonable.

There is still to do, but it has already a meaningful functionality, and performance.

Unfortunately the direct indexing of 3D arrays (that have larger overhead than 1D/2D) has a rather large overhead over the optimal strategy (10-15x) because the compiler is not smart enough to optimize it fully, but it is easy to avoid direct indexing, and as reference implementation direct indexing is still acceptable.