Segfault when calling function

[Dr Lazarus]

I receive the following error when calling the principal function in my program:

Program received signal SIGSEGV, Segmentation fault.
0x0000000000402ba8 in hoop (srccol=Cannot access memory at address 0x7fffed410cec
) at Looptrace.cc:62
62      void hoop(int srccol, int slice, int ipar, int cotable[][2][ndim], complex<double> qtrace[])

The function itself calls a number of other functions. A far as I can see, I there are not any errors with array bounds etc. The function which crashes the instant it is called is 'hoop'.Here is the program (Looptrace.cc). i have removed the last few functions because of space, but they are not used til half way through hoop, which is never reached:

#include <cstdio>
#include <cstdlib>
#include <complex>
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
using namespace std;
const int ndirac = 4;
const int ndim = 4;
const int Nc = 3;
const int Lx = 24;
const int Ly = 24;
const int Lz = 24;
const int Lt = 48;
const int spaceh = (Lx*Ly*Lz)/2;
const int volh = (Lx*Ly*Lz*Lt)/2;

void hoop(int srccol, int slice, int ipar, int cotable[][2][ndim], complex<double> qtrace[]);
void readslice(complex<double> Psi[][Nc][spaceh][ndirac], int srccol, int slice, int ipar);
complex<double> convert_complex(complex<double> in);
void gams(complex<double> a[][ndirac], complex<double> x[][ndirac], int npre, int post);
void prem(complex<double> a[][ndirac], complex<double> b[][ndirac], int npre);
void post(complex<double> x[][ndirac], complex<double> b[][ndirac], int npost);
void snum(int cotable[][2][ndim]);

int main(int argc, char *argv[])
{
  complex<double> qtrace[Lx], qsum[Lx][Lt];
  int cotable[volh][2][ndim];
  int srccol, ipar, slice, jk, jt;
  
  snum(cotable);

  for(jk = 0; jk < Lx; jk++)
  {
    for(jt = 0; jt < Lt; jt++)
    {
      qsum[jk][jt] = 0;
    }
  }
  
  for(srccol = 0; srccol < Nc; srccol++)
  {
    for(slice = 0; slice < Lt; slice++)
    {
      jt = slice;
      for(ipar = 0; ipar < 2; ipar++)
      {
	hoop(srccol, slice, ipar, cotable, qtrace);
	for(jk = 0; jk < Lx; jk++)
	{
	  qsum[jk][jt] = qsum[jk][jt] + qtrace[jk];
	}
      }
    }
  }

  return 0;
}

void hoop(int srccol, int slice, int ipar, int cotable[][2][ndim], complex<double> qtrace[])
{
  complex<double> Psi[ndirac][Nc][spaceh][ndirac];
  complex<double> qmatq[ndirac][ndirac], qmat2[ndirac][ndirac];
  complex<double> ggq[ndirac][ndirac], gg2[ndirac][ndirac];
  complex<double> qzero, qlocal, qdum;
  int source[ndim], sink[ndim], isite;
  qzero = 0;

  for(int i = 0; i < ndim; i++)
  {
    source[i] = 1; // Initial source point is always at the site (1,1,1,1)
  }

  for(int jk = 0; jk < Lx; jk++)
  {
    qtrace[jk] = qzero; // Initialize sum
  }
  
  readslice(Psi, srccol, slice, ipar);
  
  for(int jslic = 0; jslic < spaceh; jslic++)
  {
    isite = jslic*spaceh; // isite = global site number
    for(int i = 0; i < ndim; i++)
    {
      sink[i] = cotable[isite][ipar][i];
    }
    
    for(int sinkcol = 0; sinkcol < Nc; sinkcol++)
    {
      int jk;
      for(int i = 0; i < ndirac; i++) // sink spin
      {
	for(int j = 0; j < ndirac; j++) // source spin
	{
	  qmatq[i][j] = Psi[i][sinkcol][jslic][j];
	  qdum = Psi[i][sinkcol][jslic][j];
	  qmat2[j][i] = -conj(qdum);
	}
      }
    
      gams(qmat2, gg2, 5, 5);
      gams(qmatq, ggq, 3, 1);

      qlocal = qzero;
      for(int i = 1; i < ndirac; i++)
      {
	for(int j = 1; j < ndirac; j++)
	{
	  qlocal = qlocal + gg2[i][j]*ggq[j][i];
	}
      }

      jk = sink[1] + 1;
      qtrace[jk] = qtrace[jk] + qlocal;
    }
  }
  return;
}

void readslice(complex<double> Psi[][Nc][spaceh][ndirac], int srccol, int slice, int ipar)
{
  ifstream::pos_type size; // Equivalent to an integer
  complex<double> *zbuff; // Array to store binary input later in function
  int kount;

  stringstream ss;
  string filename;
  for(int srcdir = 0; srcdir < ndirac; srcdir++)
  {
    ss << "bqcd.567.00000.00." << srcdir+1 << "." << srccol << "." << slice << "." << ipar << ".prop";
    // That was the (literally) variable filename
    filename = ss.str(); // Create string
    ss.seekp(0, ios::beg); // Place the put pointer back to start of the stringstream
    cerr <<"The filename is:" << filename << endl; // A check to see which files were used
    
    ifstream readfile;
    readfile.open(filename.c_str(), ios::in | ios::binary | ios::ate);
    size = readfile.tellg(); // Determine file size fron the get pointer, which is now at the end (ate)
    readfile.seekg(0, ios::beg); // Place the get pointer at the beginning of the ifstream
    zbuff = new complex<double>[size]; // Create space for input
    readfile.read(reinterpret_cast <char *> (zbuff), size); // Read binary input into zbuff
    
    kount = 0;
    for(int k = 0; k < spaceh; k++)
    {
      for(int j = 0; j < Nc; j++) // Sink colour
      {
	for(int i = 0; i < ndirac; i++) // Sink spin
	{
	  Psi[i][j][k][srcdir] = convert_complex(zbuff[kount]); // Put complex entries into Psi
	  kount++;
	}
      }
    }

    delete zbuff;
    readfile.close();
  }
  
  cout << Psi[1][1][1][0] << endl; // Test output
  cout << Psi[2][1][2][1] << endl;
  cout << Psi[0][0][0][2] << endl;
  cout << Psi[0][0][790][3] << endl;
  
  return;
}

complex<double> convert_complex(complex<double> in) // Function to reverse the byte order of one c.number
{
  complex<double> out;
  char *p_in = (char *) &in;
  char *p_out = (char *) &out;
  
  p_out[0] = p_in[15];
  p_out[1] = p_in[14];
  p_out[2] = p_in[13];
  p_out[3] = p_in[12];
  p_out[4] = p_in[11];
  p_out[5] = p_in[10];
  p_out[6] = p_in[9];
  p_out[7] = p_in[8];
  p_out[8] = p_in[7];
  p_out[9] = p_in[6];
  p_out[10] = p_in[5];
  p_out[11] = p_in[4];
  p_out[12] = p_in[3];
  p_out[13] = p_in[2];
  p_out[14] = p_in[1];
  p_out[15] = p_in[0];
  
  return out;
}

void gams(complex<double> a[][ndirac], complex<double> x[][ndirac], int npre, int npost)
{
  complex<double> b[ndirac][ndirac];
  prem(a, b, npre);
  post(x, b, npost);

  return;

}

If anyone can see the problem I would be most grateful.

[jsmith]

You have several multi-dimensional arrays of objects on the stack in hoop(). Are you sure you aren't running out of stack space? (The fault address looks a lot like a stack address).

[Dr Lazarus]

You have several multi-dimensional arrays of objects on the stack in hoop(). Are you sure you aren't running out of stack space? (The fault address looks a lot like a stack address).

That could be the case, yes. The segmentation fault occurs the instant the function is called. It is not possible to step further into it. Also, if the issue is with large arrays, why am I receiving an error about srccol (an integer) being unaddressable? Is the stack so overloaded that the storage for all other variables has been eliminated?

[jsmith]

It could just be the debugger (gdb I assume). I get that error quite frequently when using gdb.

It would also make sense that it occurs instantly, because many of your arrays are arrays of complex<> objects that have default constructors which are going to run and attempt to initialize their data members.

[Dr Lazarus]

Yes it is gdb I'm using. Someone has suggested running the program through something called 'valgrind' although I've not heard of that before now.

It would also make sense that it occurs instantly, because many of your arrays are arrays of complex<> objects that have default constructors which are going to run and attempt to initialize their data members.

I see. The first complex array I pass (and presumably the first to be allocated on the stack) I require to be 'empty'. Does your comment mean that the default constructor initializes all its values anyway? I assume this would probably be all zeros.

Even so, that array is not particularly large (complex<double> array[24]) and is one dimensional. The big multidimensional array to be passed to the 'hoop' function is 'cotable', which is all integers. If this causes the segfault, then the memory situation is more limited than I thought. Even more disconcerting is the fact that the function 'hoop' is due to call even bigger arrays and put them on the stack, and call other functions which do the same thing. Am I fighting a hopeless battle?

[jsmith]

Just looking at your code, here are the local variables in hoop:

 complex<double> Psi[ndirac][Nc][spaceh][ndirac];
  complex<double> qmatq[ndirac][ndirac], qmat2[ndirac][ndirac];
  complex<double> ggq[ndirac][ndirac], gg2[ndirac][ndirac];
  complex<double> qzero, qlocal, qdum;
  int source[ndim], sink[ndim], isite;

ndirac = 4,
Nc = 3,
spaceh = 24 * 24 * 24 / 2 = 6912,
ndim = 4.

so Psi contains 4 * 3 * 6912 * 4 = 331,776 elements.

This probably overflows the stack.

[smilodon]

Unless you need an actual copy of the array 331,776 array, try passing by pointer or reference.

As a matter of curiosity, what does this program do?

[Dr Lazarus]

The program basically takes the trace of a matrix in the end. The end result is a number (or array of numbers) which represent something called a meson propagator in lattice QCD, a branch of physics.

By the way, I thought that arrays are passed by reference anyway? I've tested the various parts of the program, and arrays passed to each function are changed once the function is returned. Does this not mean that I wasn't simply passing a copy, but a reference?

[jsmith]

Yes. The problem is the Psi local variable inside hoop. It has 331,776 elements. Assuming an instance of complex consumes 8 bytes (2 integers * 4 bytes each for the real and imaginary parts), this data structure is over 2.5MB, which is probably going to be larger than your stack. You will need to either a) dynamically allocate the 4-D matrix or b) make it global, if possible.

a is preferred to b from an engineering standpoint; however a is not only going to be a bit of tedious coding and memory management, but also much slower than b because of all the memory allocation you need to do. This could run into performance problems if you are calling this function a lot. (I have not looked at the rest of your program).

[Dr Lazarus]

a) dynamically allocate the 4-D matrix

Presumably the stack is too limited for such a large array then. I wonder if the following code would work, or is this too naive:

complex<double> ****Psi = new complex<double> [4][3][6912][4];

I suspect this is more of a C style solution. At least it would go on the heap then, and since the machines I use have several Gb of RAM, several Mb is not a problem I would think.

[jsmith]

Yeah, that won't work. Here's a 2-d example:

int** p2dArray = new int*[ 10 ];   

for( int z = 0; z < 10; ++z )
    p2dArray[ z ] = new int[ 5 ];

p2dArray[ 8 ][ 3 ] = 4;  // eg 

[Dr Lazarus]

So for a 4-d array such as Psi[4][3][6912][4] this would be:

complex<double> ****Psi;
Psi = new complex<double>***[4];

for(int x = 0; x < 4; x++)
{
     Psi[x] = new complex<double> **[3];
     for(int y = 0; y < 3; y++)
     {
          Psi[x][y] = new complex<double> *[6912];
          for(int z = 0; z < 6912; z++)
          {
               Psi[x][y][z] = new complex<double>[4];
          }
     }
}

Am I correct?

[smilodon]

Would this work?

complex<double> *Psi = new complex<double> [4*3*6912*4];

[jsmith]

Both work. The former is most in the spirit of how people think about multi-dimensional arrays.

The latter is the method I prefer, because I inevitably end up juxtaposing indices, however
the second method requires you to write a function that converts a 4-d "coordinate" into
a linear coordinate. (Which isn't hard, and IMHO is the best solution since now the bugs
are limited to one conversion function rather than littered all over).