4
$\begingroup$

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.

$\endgroup$
  • 4
    $\begingroup$ The far larger problem here is that you are using doubly, triply, and even quadruply indirected pointers instead of a single pointer to a single memory buffer. There are no common tricks for this sort of thing, you have to do a set of loops like you're doing. This is why many scientists don't like writing code in C. $\endgroup$ – Aron Ahmadia Sep 30 '12 at 23:01
7
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ If I messed up the indexing, I will be very mad at myself. $\endgroup$ – Geoff Oxberry Oct 1 '12 at 0:05
  • $\begingroup$ I was hoping to avoid using a for loop, but I guess there's no way around it in c... $\endgroup$ – Paul Oct 1 '12 at 0:32
  • $\begingroup$ I think you could implement it using fewer for loops (again, using memcpy), but I'm not sure if you would gain in terms of performance, and you're likely to sacrifice readability and maintainability as a result. It's better to stick to for loops and worry about the details later. $\endgroup$ – Geoff Oxberry Oct 1 '12 at 0:48
  • $\begingroup$ Also, for anyone looking to this sort of example as a tutorial: I don't do any sort of error handling in the event that memory allocations fail. In production code, this possibility must be handled properly (i.e., by checking for NULL pointers returned by malloc). $\endgroup$ – Geoff Oxberry Oct 1 '12 at 2:55
5
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ This answer strikes me as more of a comment that is on the long side of being within the character limit. I also don't want to rehash a language debate. $\endgroup$ – Geoff Oxberry Oct 1 '12 at 5:45
  • $\begingroup$ That's fair. I certainly didn't intend to open the my-language-is-better-than-yours debate. It's really just a comment on the fact that some languages are better suited to some tasks than others. The question at hand here falls squarely into this category. $\endgroup$ – Wolfgang Bangerth Oct 2 '12 at 4:43
  • $\begingroup$ I agree with Wolfgang. I need this kind of thing very often, so it's good to point out to other languages that can get the job done easily. $\endgroup$ – Ondřej Čertík Oct 3 '12 at 4:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.