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):

#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;

}

void prem(complex<double> a[][ndirac], complex<double> b[][ndirac], int npre)
{
  complex<double> qI(0.0, 1.0);
  
  for(int i = 0; i < ndirac; i++)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[i][j] = 0;
    }
  }

  if(npre == 1)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[0][j] =  qI*a[3][j];
      b[1][j] =  qI*a[2][j];
      b[2][j] = -qI*a[1][j];
      b[3][j] = -qI*a[0][j];
    }
  }

  if(npre == 2)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[0][j] =  a[3][j];
      b[1][j] = -a[2][j];
      b[2][j] = -a[1][j];
      b[3][j] =  a[0][j];
    }
  }

  if(npre == 3)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[0][j] =  qI*a[2][j];
      b[1][j] = -qI*a[3][j];
      b[2][j] = -qI*a[0][j];
      b[3][j] =  qI*a[1][j];
    }
  }

  if(npre == 4)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[0][j] =  a[0][j];
      b[1][j] =  a[1][j];
      b[2][j] = -a[2][j];
      b[3][j] = -a[3][j];
    }
  }

  if(npre == 5)
  {
    for(int j = 0; j < ndirac; j++)
    {
      b[0][j] =  a[2][j];
      b[1][j] =  a[3][j];
      b[2][j] =  a[0][j];
      b[3][j] =  a[1][j];
    }
  }

  return;
}

void post(complex<double> x[][ndirac], complex<double> b[][ndirac], int npost)
{
  complex<double> qI(0.0, 1.0);
  
  for(int i = 0; i < ndirac; i++)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[i][j] = 0;
    }
  }

  if(npost == 1)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[0][j] = -qI*b[3][j];
      x[1][j] = -qI*b[2][j];
      x[2][j] =  qI*b[1][j];
      x[3][j] =  qI*b[0][j];
    }
  }

  if(npost == 2)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[0][j] =  b[3][j];
      x[1][j] = -b[2][j];
      x[2][j] = -b[1][j];
      x[3][j] =  b[0][j];
    }
  }

  if(npost == 3)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[0][j] = -qI*b[2][j];
      x[1][j] =  qI*b[3][j];
      x[2][j] =  qI*b[0][j];
      x[3][j] = -qI*b[1][j];
    }
  }

  if(npost == 4)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[0][j] =  b[0][j];
      x[1][j] =  b[1][j];
      x[2][j] = -b[2][j];
      x[3][j] = -b[3][j];
    }
  }

  if(npost == 5)
  {
    for(int j = 0; j < ndirac; j++)
    {
      x[0][j] =  b[2][j];
      x[1][j] =  b[3][j];
      x[2][j] =  b[0][j];
      x[3][j] =  b[1][j];
    }
  }

  return;
}

void snum(int cotable[][2][ndim])
{
  int ix, iy, iz, it, ipar, neven, nodd, isite;
  int site;

  neven = 0;
  nodd = 0;

  for(it = 0; it < Lt; it++)
  {
    for(iz = 0; iz < Lz; iz++)
    {
      for(iy = 0; iy < Ly; iy++)
      {
	for(ix = 0; ix < Lx; ix++)
	{
	  ipar = (ix + iy + iz + it) % 2;
	  
	  if(ipar == 0)
	  {
	    neven++;
	    isite = neven;
	  }
	  else if(ipar == 1)
	  {
	    nodd++;
	    isite = nodd;
	  }
	  
	  cotable[isite][ipar][0] = ix;
	  cotable[isite][ipar][1] = iy;
	  cotable[isite][ipar][2] = iz;
	  cotable[isite][ipar][3] = it;
	  
	}
      }
    }
  }
  return;
}

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

[Lindley]

You don't have any check in snum() that isite < volh. Are you certain that this will always be the case? An assert() might be appropriate.

Other than that, run it through valgrind if you're on Linux. That will pinpoint when things start to go wrong.

[Dr. Lazarus]

Thanks for your reply. I used asserts as you suggested, and they do in fact fail. I've just used valgrind with the following message:

 Warning: client switching stacks?  SP change: 0x7FF0003A0 --> 0x7FE5DB9B0
==21398==          to suppress, use: --max-stackframe=10635760 or greater
==21398== Invalid write of size 4
==21398==    at 0x40340F: main (Looptrace.cc:29)
==21398==  Address 0x7fe5db9fc is on thread 1's stack
==21398==
==21398== Process terminating with default action of signal 11 (SIGSEGV)
==21398==  Access not within mapped region at address 0x7FE5DB9FC
==21398==    at 0x40340F: main (Looptrace.cc:29)
==21398==
==21398== Invalid write of size 8
==21398==    at 0x4802338: _vgnU_freeres (vg_preloaded.c:56)
==21398==  Address 0x7fe5db9a8 is on thread 1's stack
==21398==
==21398== Process terminating with default action of signal 11 (SIGSEGV)
==21398==  Access not within mapped region at address 0x7FE5DB9A8
==21398==    at 0x4802338: _vgnU_freeres (vg_preloaded.c:56)
==21398==
==21398== ERROR SUMMARY: 2 errors from 2 contexts (suppressed: 5 from 1)

Am I right in saying that an invalid write of size 4 and 8 is for an int and a double? Looptrace.cc:29 is the first line of main. In fact it's the first bracket. By the way, why is valgrind telling me that the memory problem is right at the start of main, when it occurs at the call of hoop? Using the --num-callers option to show more of the stack trace has no effect on this message.

[Lindley]

Did you compile with the -g flag (debug symbols) and *without* any optimization? That could affect the line numbers it gives.

By the way, why is valgrind telling me that the memory problem is right at the start of main, when it occurs at the call of hoop?

Often the cause of a memory-related problem will occur well before the program actually implodes. The segfault occurs when things finish collapsing; valgrind is able to give you a much better idea of when they start, which is usually much more useful for solving the problem.

[Dr. Lazarus]

Did you compile with the -g flag (debug symbols) and *without* any optimization? That could affect the line numbers it gives.

Yes, I compiled with -g and no optimization. By the way, someone has suggested I should pass some of these arrays by reference, but I had the impression that this is what happens anyway. These arrays behave as though they were passed by reference, since they are modified outside the relevant functions.

[Lindley]

What's actually happening is that you're passing a pointer to the start of the array by value. Arrays in C/C++ aren't objects, just memory ranges, so technically it's impossible to pass them at all, only pointers to them. (std::vectors are objects though.)

I thought of a possible cause. You're creating a lot of arrays on the stack at the start of main, several of them multi-dimensional. In particular, cotable will be very large! (Over 10 megabytes.)

The stack isn't all that large in C++, maybe a couple of megabytes. You shouldn't construct really large arrays on the stack----you may simply be overflowing it. Either create cotable on the heap using new[] (and don't forget to delete[] it!), or change it into a std::vector.

[Dr. Lazarus]

I might do this with raw pointers on the heap, even though that will be fiddly. I think you are right, I'm asking (far) too much of the stack.