This HTML automatically generated with rman for NEMO
Table of Contents


mdarray - hierarchical multi-dimensional (C style row-major) array allocators with sequential memory


#include <mdarray.h>
#define MDMAXDIM    8
typedef real     *mdarray1;   /* vector */typedef
mdarray1 *mdarray2;   /* matrix */typedef mdarray2 *mdarray3;   /* tensor
*/typedef mdarray3 *mdarray4;   /* hyper tensor */...typedef mdarray7 *mdarray8;
  /* see MDMAXDIM */
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, i.e. row-major. Up to MDMAXDIM (defined via the header file) dimensions are available. Actual data memory is guarenteed to be sequentually in memory, as they are with static multidimensional arrays. This makes it easy to overlay them with traditional languages such a Fortran (with the caution of row-major and column-major differences).

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, but opposite those of fortran. The matching fortran array would be A(10,20) in this case.

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)

When using this with 3rd party software that needs to write into contiguous memory, the address of the first element will normally be sufficient, e.g.

      mdarray2 data2 = allocate_mdarray2(nrows,ncols);
      fits_read_col(fptr, TDOUBLE, data_col, 1, 1, ncols, &nulval, &data2[0][0],
&anynul, &status);      


The TESTBED function (mdarraytest) exercizes 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

vectmath(3NEMO), mdbench(1NEMO)
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
dec-2019    Increased MDMAXDIM to 8        PJT

Table of Contents