Extracting a 2D matrix from a multidimensional array in C

Question

In my c language program, I have to store multiple dense $m\times m$ matrices corresponding to gridpoints $x_i$ with $i=1,...,n$. I decided to create a three dimensional array $A\in R^{n\times m\times m}$ as follows:

double ****A;
A=(double****)malloc(sizeof(double***)*n);  
for (i=0; i<n; i++) {
     A[i]=(double***)malloc(sizeof(double**)*m);
     for (j=0; j<m; j++) {
          A[i][j]=(double**)malloc(sizeof(double*)*m);
     }
}

In this way, I'm indexing my matrix as A[i][j][k] such that at each index [i] corresponds to a unique $m\times m$ matrix. However, I now have to input these 2 dimensional ($m\times m$)matrices into a solver which assumes they are stored as only two dimensional arrays.

If i were writing the program in matlab, I could extract the 2D matrix by simply coding:

B=A(i,:,:);

But I'm not sure how to do this in c. There must be some clever trick using pointers.

Answer 1

As Aron notes in his comment, rarely are matrices stored in nested pointer data structures in practice. Multiple levels of indirection (I'm told) cause considerable performance penalties, and require multiple allocations and frees (more performance issues, plus more potential to screw up and segfault). I've seen nested pointers used to store matrices only a few times (in CHEMKIN, and I think also in Cantera). Far more common practice is to store the data in memory associated with a single pointer; the C interfaces to BLAS and LAPACK use this practice.

You probably want to do something more like

double *A;
A = (double *) malloc(sizeof(double) * m * m * n);
// A[(i - 1) * m * m + (j - 1) * m + k] == A[i][j][k]

where i, j, k refer to traditional one-based mathematical indices; this sort of arrangement of data is called row-major ordering.

Then, to do the slicing you'd like, you will need to copy the data in A to another pointer:

double *B;
B = (double *) malloc(sizeof(double) * m * m);
// Here, zero-based indexing is used instead of the one-based indexing in the
// last code snippet.
for (j = 0; j++; j < m) {
    for (k =0 ; k++; k < m) {
        B[j * m + k] = A[i * m * m + j * m + k]; // B == A[i, :, :]
    }
}

There might be slicker ways of doing this operation, using something like memcpy, but you'll still need to copy the data to another data structure.

Let's suppose that you proceed with your current code, and suppose further that you continue with your strategy of storing a matrix in nested pointers. To slice A, you would do something like:

double **B;
B = (double **) malloc(sizeof(double *) * m);
// Zero-based indexing also used in this code snippet.
for (j = 0; j++; j < m) {
    B[j] = (double *) malloc(sizeof(double *) * m);
    for (k = 0; k++; k < m) {
        B[j][k] = A[i][j][k];
    }
}

Again, there might be slicker ways of doing this copying, such as using memcpy.

In addition to MATLAB, languages such as Fortran 90/95/2003/2008 and Python support array slicing syntax natively.

Answer 2

If you don't want to litter your program with this kind of copying loop then maybe C isn't your language. Assuming this is a common operation, C++ provides far nicer ways to do this sort of thing (by making the array a class and, possibly, making the sub-array a view of it). Even Fortran has a nicer syntax for this kind of operation. What this illustrates is that you can express everything in every Turing-complete language but that it isn't true that everything can be expressed with equal ease in each language.