How to allocate array starting negative index

Question

I am trying to allocate a 3D array u[-nx/2:nx/2-1][-nx/2:nx/2-1][-nx/2:nx/2-1]

    int nx = 512;
    double *** u = (double ***)malloc(nx * sizeof(double**));

    for (int i = -nx/2; i < nx/2; i++) {
        u[i] = (double **)malloc(nx * sizeof(double *));
        for (int j = -nx/2; j < nx/2; j++) {
            u[i][j] = (double *)malloc(nx * sizeof(double));
        }
    }

Is this a correct way to do it? If it's not, how should I change it?

Answer 1

First off, all C arrays have index values ranging from 0 to ELEMENT_COUNT-1. Always.

As you are using malloc, I am presuming that the value of nx is only defined at runtime. This rules out static array sizes and thus rules out using the cute arr[x][y][z] syntax as in:

#define NX 512

double arr[NX][NX][NX];

void foo(void)
{
   ...
   arr[z1][y1][x1] += 2 * arr[z2][y2][x2];
   ...
}

That in turn means that to have the functionality of a 3D array with nx different values for each of its three dimensions dimension, you will need to allocate a linear memory area of size nx_cubed = nx * nx * nx. To calculate that value nx_cubed properly, you will need to check for integer overflows.

Also, you need to properly convert from signed int coordinate values to unsigned size_t values used in the 0 based index ranges.

if (nx < 0) {
    fprintf(stderr, "negative value of nx\n");
    exit(EXIT_FAILURE);
}

const size_t unx = nx;
const size_t nx_cubed = unx * unx * unx;
/* TODO: Complete check for overflows */
if (nx_cubed < unx) {
    fprintf(stderr, "nx_cubed overflow\n");
    exit(EXIT_FAILURE);
}

Then you can allocate a memory buffer of the appropriate size, and then check that the malloc call has actually worked.

double *buf = malloc(nx_cubed);
if (!buf) {
    fprintf(stderr, "Error allocating memory for nx_cubed elements\n");
    exit(EXIT_FAILURE);
}

Now there is the question of calculcating the array index from your x, y, and z values each ranging from -nx/2 to nx/2-1. I recommend writing a function for that which maps that range to the 0 to nx-1 range, and then calculates the proper linear index from the three 0-based values. Again, proper integer overflow checks should be performed.

size_t array3index(const size_t nx, const int x, const int y, const int z) {
   const size_t half_nx = nx/2;
   /* zero based 3D coordinates,
    * this probably triggers some signedness warnings */
   const size_t x0 = half_nx + x;
   const size_t y0 = half_nx + y;
   const size_t z0 = half_nx + z;
   if ((x0 >= nx) || (y0 >= nx) || (z0 >= nx)) {
     fprintf(stderr, "Signed coordinate(s) out of range: (%d, %d, %d)\n",
             x, y, z);
     exit(EXIT_FAILURE);
   }
   const size_t idx = nx * (nx * z0 + y0) + x0;
   /* Assuming that we have already checked that nx*nx*nx does not
    * overflow, and given that we have checked for x0, y0, z0 to be
    * in the range of 0 to (nx-1), the idx calculation should not
    * have overflown here. */
   return idx;
}

Then you can do your accesses to the 3D array like

const i1 = array3index(nx, x1, y1, z1);
const i2 = array3index(nx, x2, y2, z2);

buf[i1] += 2*buf[i2];

Considering the amount of calculations needed inside array3index, I would examine whether it makes more sense to do the array iteration in the 0 to nx-1 domain directly, and only convert that to -nx/2 to nx/2-1 range values if you actually need that value within a calculation.

Answer 2

No. Array in C is actually just plain/flat memory block, which is always 0 based and always in 1d (one demension).

Suppose you need a 3d array in arbitrary boundary, say u[lb_1d, ub_1d][lb_2d, ub_2d][lb_3d, ub_3d], you will need to do some mapping -- address space from 3d to 1d and vice versa --.

Sample implementation like this:

typedef struct
{
    double* _arr;
    int _lb_1d;
    int _ub_1d;
    int _lb_2d;
    int _ub_2d; 
    int _lb_3d;
    int _ub_3d;
}DoubleArr3D;

DoubleArr3D*  create_3d_arr(int lb_1d, int ub_1d, int lb_2d, int ub_2d, int lb_3d, int ub_3d) 
{
    int array_size = (ub_1d - lb_1d +1) * (ub_2d - lb_2d +1) * (ub_3d - lb_3d +1);
    DoubleArr3D * arr = (DoubleArr3D *)malloc( sizeof( DoubleArr3D) );
    if (!arr)
    {
        return NULL;
    }

    arr->_lb_1d = lb_1d;
    arr->_ub_1d = ub_1d;
    arr->_lb_2d = lb_2d;
    arr->_ub_2d = ub_2d;
    arr->_lb_3d = lb_3d;
    arr->_ub_3d = ub_3d;

    arr->_arr = (double*) malloc(sizeof(double) * (size_t) array_size);
    if (!arr)
    {
        free(arr);
        return NULL;
    }

    return arr;
}

// arr[i1d, i2d, i3d] ==> arr_get_at(arr, i1d, i2d, i3d)
double arr_get_at(DoubleArr3D*  arr, int i1d, int i2d, int i3d )
{
    if (!arr || !arr->_arr)
    {
        // just demo of 'validation check'. in real code we should have meanful error report
        return 0;
    }

    return arr->_arr [  
            i3d - arr->_lb_3d 
            + (i2d - arr->_lb_2d ) * (arr->_ub_3d - arr->_lb_3d +1) 
            + (i1d - arr->_lb_1d ) * (arr->_ub_2d - arr->_lb_2d +1) * (arr->_ub_3d - arr->_lb_3d +1)
    ];
}

Answer 3

No, that’s not correct. You can get it to work by placing every pointer in the middle of the dimension it represents:

int nx = 512;
double*** u = (double***)malloc(nx * sizeof(double**)) + nx/2;

for (int i = -nx/2; i < nx/2; i++) {
    u[i] = (double**)malloc(nx * sizeof(double*)) + nx/2;
    for (int j = -nx/2; j < nx/2; j++) {
        u[i][j] = (double*)malloc(nx * sizeof(double)) + nx/2;
    }
}

but that’s unusual and confusing, does a lot of separate allocations, and has to be undone for the deallocation step.

Consider one block with accessors instead:

#define NX 512

/* or just double* if nx is dynamic, and calculate the index manually */
double[NX][NX][NX]* u = malloc(sizeof(*u));

double array_get(double[NX][NX][NX] const* u, int i, int j, int k) {
    return u[i + NX/2][j + NX/2][k + NX/2];
}

void array_set(double[NX][NX][NX]* u, int i, int j, int k, double value) {
    u[i + NX/2][j + NX/2][k + NX/2] = value;
}