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

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

Code:
#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.