#include <mdarray.h> #define MDMAXDIM 8 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);...
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 fortan. 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);
% 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+0wOn grus (P3/930) 5.0 vs. 19.9 sec. On chand (P4/2500) 1.1 vs. 12.3 sec CPU.
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)
vectmath(3NEMO), mdbench(1NEMO) nrutil.c (Numerical Recipes) for an alternative approach GSL: http://www.gnu.org/software/gsl/manual/html_node/Matrices.html (cf. valarray in C++) http://www.fredosaurus.com/notes-cpp/arrayptr/23two-dim-array-memory-layout.html http://www.fftw.org/ http://c-faq.com/aryptr/
~/src/kernel/misc mdarray.c
5-may-03 Created PJT 19-feb-06 Allocate actually memory sequentually PJT dec-2019 Increased MDMAXDIM to 8 PJT