Table of Contents


mdarray - simple multi-dimensional (C style) array allocators with sequential memory


#include <mdarray.h>
#define MDMAXDIM    7 
typedef real     *mdarray1;   /* vector */typedef
mdarray1 *mdarray2;   /* matrix */typedef mdarray2 *mdarray3;   /* tensor
*/typedef mdarray3 *mdarray4;   /* hyper tensor */...
mdarray1 allocate_mdarray1(n1)mdarray2
allocate_mdarray2(n2,n1)mdarray3 allocate_mdarray3(n3,n2,n1)mdarray4 allocate_mdarray4(n4,n3,n2,n1)...void
free_mdarray1(mdarray1 x, int n1);void free_mdarray2(mdarray2 x, int n2,
int n1);void free_mdarray3(mdarray3 x, int n3, int n2, int n1);void free_mdarray4(mdarray4
x, int n4, int n3, int n2, int n1);...


mdarray offers a uniform way to allocate and free multi-dimensional (real) C-style arrays. Up to MDMAXDIM (defined via the header file) dimensions are available. Actually data memory is guarenteed to be sequentually in memory, as they are with static multidimensional arrays.

Arrays are stored in row major order, if rows vary most rapidly in memory. For example in

   double matrix a[20][10];
the right-most index (10) varies most rapidly in memory, exactly like static C arrays.

Example of use: ( b = A.x)

    mdarray1 x = allocate_mdarray1(10);      /*  x[10]     */
    mdarray2 A = allocate_mdarray2(20,10);   /*  A[20][10] */
    mdarray1 b = allocate_mdarray1(20);      /*  b[20]     */
    int i,j;
    for (j=0; j<20; j++) {
        b[j] = 0.0;
        for (i=0; i<10; i++)
            b[j] += A[j][i]*x[i]

When referring to rows and columns, we write these as A[row][col] (e.g. in C, or the current mdarray2) or A(row,col) (e.g. in Fortran)


The TESTBED function (mdarraytest) excersizes a simple initialization of all array elements, and summing them up. This also clearly shows how reversing the order of array access affects the CPU speed , viz. timings on a P4/1600 laptop:
    % time mdarraytest 4,4,4 iter=1000000  flip=f
    Working with ndim=3 MD_Array[[4][4][4]
    2.110u 0.010s 0:02.31 91.7%     0+0k 0+0io 135pf+0w
    % time mdarraytest 4,4,4 iter=1000000  flip=t
    Working with ndim=3 MD_Array[[4][4][4]
    2.050u 0.000s 0:02.14 95.7%     0+0k 0+0io 135pf+0w
    % time mdarraytest 100,100,100 iter=100 flip=f
    Working with ndim=3 MD_Array[[100][100][100]
    3.210u 0.020s 0:03.35 96.4%     0+0k 0+0io 135pf+0w
    % time mdarraytest 100,100,100 iter=100 flip=t
    Working with ndim=3 MD_Array[[100][100][100]
    30.530u 0.040s 0:31.85 95.9%    0+0k 0+0io 135pf+0w
On grus (P3/930) 5.0 vs. 19.9 sec. On chand (P4/2500) 1.1 vs. 12.3 sec CPU.


Although both static and dynamic multi-dimensional arrays
    real     s3[4][5][6];
    mdarray3 d3 = allocate_mdarray3(4,5,6);
use sequential memory, the dynamic array uses extra memory to use the arrays of arrays, and therefore the addresses
    d3 !=  d3[0] != d3[0][0]
are not the same, whereas for static arrays they are
    s3 ==  s3[0] == s3[0][0]
For example, for a double precision a[n3][n2]n1] mdarray3, instead of using 2*n1*n2*n3 words, it will use ((2*n1+1)*n2+1)*n3 words.

To quote the FFTW manual: using it will cause FFTW to die a painful death (referring to such dynamic arrays)

See Also

nrutil.c (Numerical Recipes) for an alternative approach
GSL: (cf. valarray
in C++)


Peter Teuben


~/src/kernel/misc    mdarray.c

Update History

5-may-03    Created       PJT
19-feb-06    Allocate actually memory sequentually    PJT

Table of Contents