Bus error using mpic++ but not using g++

Question

I have some code written in C++ that compiles and runs fine (with reasonable output) when using g++, but when I try to use mpic++, I get a runtime bus error. I've managed to figure out where the bus error is occuring, but not why. Here is my code:

one = (double *) malloc(sizeof(&one) * nx * nx * nt);
two = (double **) malloc(sizeof(&two) * nx * nx);
rho_exp = (double ***) malloc(sizeof(&rho_exp) * nx);

for(i = 0; i < nx * nx; i++)
    two[i] = &one[i * nt];

for(i = 0; i < nx; i++)
    rho_exp[i] = &two[i * nx];

for(i = 0; i < nx; i++)
    for(j = 0; j < nx; j++)
        for(k = 0; k < nt; k++)
            rho_exp[i][j][k] = 0;

The bus error is occurring during the three nested for loops. My question is two-fold: One, I assume, I've screwed up my allocation for the 3D matrix; how should I have done it? Two, why did this work with gcc and g++, but not with mpic++?

Answer 1

One thing that jumps out is that you are using sizeof(&foo) (size of a pointer to foo) when you probably mean sizeof(*foo) (size of the thing that foo points to).

I bet what is happening is that mpic++ is compiling to a 32-bit target, where the size of doubles is 64 bits and the size of addresses is 32 bits. That gives you a mismatch and causes a problem. Your g++ is probably targeting a 64-bit system, where those sizes are the same (both 64 bits).

EDIT:

Here's the code you should be using:

  double * one = new double[nx * nx * nt];
  double ** two = new double*[nx * nx];
  double ***rho_exp = new double**[nx];

And in all honesty, you are reinventing the wheel. You should be using your vendor-provided BLAS library, which has fast matrix operations optimized for your platform (probably hundreds of times faster than whatever you write...).

Answer 2

Here's the general procedure for dynamically allocating an X x Y x Z array:

double ***rho_exp = malloc(sizeof *rho_exp * X);
if (rho_exp)
{
  size_t i;
  for (i = 0; i < X; i++)
  {
    rho_exp[i] = malloc(sizeof *rho_exp[i] * Y);
    if (rho_exp[i])
    {
      size_t j;
      for (j = 0; j < Y; j++)
      {
        rho_exp[i][j] = malloc(sizeof *rho_exp[i][j] * Z);
        if (rho_exp[i][j])
        {
          size_t k;
          for (k = 0; k < Z; k++)
          {
            rho_exp[i][j][k] = 0;
          }
        }
      }
    }
  }
}

Never mind, that's for C; since you're using C++, use new as described above.