Heap exhausted far too quickly

[Dr Lazarus]

Some may recognize this program from earlier posts - my principal problem is to allocate very large arrays for use in the program. The two large arrays are int cotable[331,776][2][4] and complex<double> Psi[4][3][6912][4], which is 2,654,208 and 331,776 elements respectively. After these are (successfully) allocated, I do some binary file input which requires me to allocate a char array (called buffer). Cotable is ~10.6 Mb, Psi is ~5.3 Mb, buffer is ~1.3 Mb.

My problem is this: The use of 'new' to allocate cotable and Psi is successful. However, when the program reaches buffer = new char[size] I receive a bad_alloc. How can the heap be so limited? My machine has several Gb of system memory and, presumably, even more virtual memory. The sum total of these arrays is not more than 18 Mb. Can anyone see what the problem is? I have removed the last few functions due to space limitations:

#include <cstdio>
#include <cstdlib>
#include <complex>
#include <math.h>
#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <new>
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, complex<double> ****Psi, complex<double> qtrace[]);
void readslice(complex<double> ****Psi, 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);
void out_of_store();

int main(int argc, char *argv[])
{
  complex<double> qtrace[Lx], qsum[Lx][Lt];
  int srccol, ipar, slice, jk, jt;
  
  set_new_handler(out_of_store);

  // The following code allocates memory on the heap for cotable[volh][2][ndim], a very large array
  int ***cotable; //Pointer to 2-d arrays (which are pointers to pointers **)
  cotable = new int **[volh]; // Cotable now points to volh different two-d arrays. It is a 3-d array.

  if(cotable == NULL)
  {
    //Balls! This pointer points into outer space, where no pointer has gone before
    cout << "cotable pointer is null. Oops!" << endl;
  }
  
  for(int x = 0; x < volh; x++)  //Assigning memory for the other two dimensions
  {
    cotable[x] = new int *[2];
    for(int y = 0; y < 2; y++)
    {
      cotable[x][y] = new int[ndim];
    }
  }

// The following code allocates memory for Psi. Move to main if there's a problem
  complex<double> ****Psi;
  Psi = new complex<double> ***[ndirac];

  if(Psi == NULL)
  {
    //Yikes! Pointer not pointing anywhere...
    cout << "Psi pointer is null. Oops!" << endl;
  }

  for(int x = 0; x < ndirac; x++)
  {
    Psi[x] = new complex<double> **[Nc];
    for(int y = 0; y < Nc; y++)
    {
      Psi[x][y] = new complex<double> *[spaceh];
      for(int z = 0; z < spaceh; z++)
      {
	Psi[x][y][z] = new complex<double> [ndirac];
      }
    }
  }
  
  
  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++)
    {
      for(ipar = 0; ipar < 2; ipar++)
      {
	hoop(srccol, slice, ipar, cotable, Psi, qtrace);
	for(jk = 0; jk < Lx; jk++)
	{
	  qsum[jk][slice] = qsum[jk][slice] + qtrace[jk];
	}
      }
    }
  }

  return 0;
}

void hoop(int srccol, int slice, int ipar, int ***cotable, complex<double> ****Psi, complex<double> qtrace[])
{
  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];
    }
    
    // Code to delete cotable
    for(int x = 0; x < volh; x++)
    {
      for(int y = 0; y < 2; y++)
      {
	delete[] cotable[x][y];
      }
    }
    for(int x = 0; x < volh; x++)
    {
      delete[] cotable[x];
    }
    delete[] cotable;

    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, int srccol, int slice, int ipar)
{
  //ifstream::pos_type size; // Equivalent to an integer
  long size; // Integer to record file size
  complex<double> *zbuff; // Pointer to complex to create complex array
  //char *buffer; // Pointer to char - required to store binary data, which is a whole bunch of chars
  int kount;

  stringstream ss;
  string filename;
  for(int srcdir = 0; srcdir < ndirac; srcdir++)
  {
    ss << "bqcd.567.00000.00." << srcdir+1 << "." << srccol+1 << "." << 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 from 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
    //buffer = new char[size]; // Create space for input
    char buffer[size];
    readfile.read(buffer, size); // Read binary input into buffer, which is an array of chars
    zbuff = reinterpret_cast <complex<double> *> (buffer); // Store the base address of buffer
    // in a pointer to complex numbers. The data is now an array of complex numbers
    
    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++;
	}
      }
    }

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

[jsmith]

I don't see a new char[] anywhere in your program, so I can't be sure.

Make sure that "size" is a reasonable value. Also are you calling these functions several times? One possibility is that you are running out of virtual memory (NOT swap space) in your process' address space due to fragmentation. Arrays, when allocated, by definition are contiguous in memory. If your process does not have a contiguous block of memory large enough for the allocation, bad_alloc would be thrown. The easiest way to determine this is to print out the pointer returned by new everywhere you use new. Fragmentation will result in pointer values increasing over time toward the top of the process' address space.

Two other things:
1) Unless you are using the nothrow version of new (you aren't), there is no point to checking the returned pointer against NULL; the throwing version of new will never return NULL.

2) What you are doing in convert_complex is extremely dangerous. You are effectively memcpy()ing an object (except in your case, it is more like a reverse memcpy).

[Dr Lazarus]

Apologies, the new char[size], just commented out, it is at line 205. This line is called multiple times, but when I call new char[] I also delete 'buffer' before the end of the loop in 'readslice'.

Unfortunately I was forced to write the function convert_complex since the binary files I read were not created by me, and there is no suitable gcc to reverse the endianness on my system. The function does work, since the numbers I acquire now make sense.

The program as it stands with the commented out parts in readslice returns a segfault, but really I want it with new char[size] called, and deleted, each time in the loop, but this is the line which returns bad_alloc. It doesn't make any sense that the this would be a problem, since I've already allocated much larger items earlier in the program.

[helios]

Did you check that all the expressions inside the [] for new are greater than zero? Between lines 203 and 205, for example, there are no checks that the file is not empty.

[Dr Lazarus]

I seem to have identified the problem and you are correct, there is a trivial error in my construction of the filename - I need to pad zeroes when I insert my variable 'slice'. If I don't do this, the file does not open, and size = 0; I'm surprised that the bad_alloc complaint does not happen at that point in the code, but such weirdness is common I believe.