Dynamic Contiguous 3D arrays in C

Question

I'm trying to implement dynamically allocated contiguous 3D arrays in a C code. The arrays must be contiguous because I'm relying on netCDF output of the arrays. Now I adapted a solution posted here Stack OverFlow Solution. This works fine for dynamically allocating the arrays and indexing them...however, when netCDF outputs them there is an offset which appears to scale as the 2nd index size (jcount). Here is the modified function:

void*** newarray(int icount, int jcount, int kcount, int type_size)
{
  int i,j,k;
  void*** iret = (void***)malloc(icount*sizeof(void***)+icount*jcount*sizeof(void**)+icount*jcount*kcount*type_size);
  void** jret = (void**)(iret+icount);
  char* kret = (char*)(jret+icount*jcount);
  for(i=0;i<icount;i++)
  {
    iret[i] = &jret[i*jcount];
  }
  for(i=0;i<icount;i++)
  {
    for(j=0;j<jcount;j++)
    {
      jret[i*jcount+j] = &kret[i*jcount*kcount*type_size+j*kcount*type_size];
    }
  }
  return iret;
}

If I understand this function correctly iret is allocated space for the 3D pointers composing iret (first index), the 2D pointers composing jret (2nd index) and the space for the actual values composing kret. The 2D jret pointer is then associated with the 2D array section of iret. Then the same is done for kret. Then every address of iret is pointed to the first value of each section of jret. Then each address of jret is pointed to the first address of kret.

For the record, everything was working fine with preprocessor defined values for my arrays. Also, if I use some printf statements in the code to check the numerics of the arrays they all appear to index correctly and the code runs correctly, it's just the output seems to be the result of non-contiguous memory storage of the arrays.

I have a structure of the form:

typedef struct
{
  double ***test;
} STRUCT_TYPE;

Which I then allocate using

mhd    = (STRUCT_TYPE *)  malloc(sizeof(STRUCT_TYPE));
mhd.test = (double***) newarray(101,7,101,sizeof(double));

This may be an issue with netCDF...but I'd just like to know my allocation routine isn't the issue.

Answer 1

Well I found a fix...although it's more of a convention thing. When I originally designed my code I simply used a preprocessor directive to declare my array sizes. Then I simply declared each array inside the structure. In order to save memory I simply passed a pointer to the structures to my sub-functions. This meant if I was going to reference an element outside of the main() routine I used something like:

grid->x;

Now and I references each element of x like

grid->x[i][j][k];

Of course, when a netCDF function needed a pointer to the variable it was going to store, I just passed it

grid->x

However, using this new memory allocation function, when I pass this to a function I'm actually passing the first address of the memory space set aside for the various pointers. In essence I was passing a pointer to the wrong address in memory. All I had to do is change the passed argument to:

&grid->x[0][0][0]

and I get the proper array output. This makes me think that when C allocates memory for an array it allocates the memory used to store the values first then it allocates the memory for the pointers which provide structure to the multidimensional array.

Maybe someone can clear this one up for me?

Answer 2

This is actually an answer to the question in Lazer's answer:

This makes me think that when C allocates memory for an array it allocates the memory used > to store the values first then it allocates the memory for the pointers which provide structure to the multidimensional array.

Maybe someone can clear this one up for me?

There is a very fundamental difference in C between

double arr[A][B][C];

and

double ***ptr;

Although both can be indexed with x[i][j][k], the resulting machine code is very different.

double arr[A][B][C] is an array of A arrays of B arrays of C doubles. When indexing this in the code as arr[i][j][k], the compiler translates that to an offset of ((i*B + j)*C + k)*sizeof(double) bytes from the start of arr. There are no intermediate pointers involved here, but the dimensions B and C must be known to the compiler.

'double ***ptr' is a pointer to (the start of an array of) pointers to (the start of an array of) pointers to (the start of an array of) double. When indexing this in the code as ptr[i][j][k], the compiler has no choice but to do each indexing operation separately and follow all the intermediate pointers. This results in code similar to

temp1 = ptr[i];
temp2 = temp1[j];
result = temp2[k];

The newarray function carves all these component-arrays out of a single block of memory, but that is not required by the language and the compiler does not have the knowledge that this is taking place, so it has to treat all the component-arrays as being completely independent of each other.