Run out of memory!

[zaratustra]

Hello, I was recently running a self-made programm in C++ and I got the following message during execution in the command window:

terminate called after throwing an instance of 'std::bad_alloc'
what(): std::bad_alloc

Afterwards I receive a Windows message saying that the programm stopped because I have not enough memory. Thing is I am using dynamic memory for the program and I am afraid I have done something wrong with it... Is there any change to fix this? thanks in advance

[Chervil]

You'd need to show the code, to see what might be the cause. If the program is very long, show a minimal version with enough code so it can be compiled and demonstrate the problem.

[zaratustra]

Ok the thing is I am using a header file for a matrix class, so the program consists on 3 files: header, functions of header and main. The program works with different loops and now, when I try to run it even for small matrix (1000 elements), it stops working, I think because the programm cannot allocate memory for the matrices anymore... Does that mean I have lost all of my RAM memory and there is no way I can get it back? I start to think this sounds pretty serious if that is the case

I will put the main part of the code now separately:
main:

  #include <iostream>
#include <cmath>
#include <fstream>
#include "matrix.h"
#include <stdlib.h>
#include <time.h>//Header of time to initiate the seed
#define pi 3.14159265359897

using namespace std;

// The main functions are defined at the beginning of the 
//program. The header file "matrix.h" includes the class matrix to work with matrices by using different operations


double* v_rand_MB(int npart, double a){ //arguments are # of particles and parameter "a", related to temperature of the gas
 int count;//counts the total number of accepted random numbers from the distribution
 count=0;
//Pointer to the array of velocities
 double* myp;
 double interval,bb,vmax,Pd,Ptest,num,den,alp;
 
 vmax=5.0*a;
 
 interval=abs(vmax);
 
 
 myp=new double[npart];//The pointer points towards an npart dimensional array of type double
  
  //loop for generation of random numbers
  while(count<npart){
  	 
  	 //Generate a random number on the desired interval.
  	 bb=rand()/(1.0*RAND_MAX);
  	 bb=interval*bb;
    	 
  	 //Probability density, i.e., MB distribution
  	 den=sqrt(2.0/pi)/(1.0*pow(a,3));
  	 alp=1.0/(2*pow(a,2));
  	 num=pow(bb,2)*(exp(-alp*pow(bb,2)));
  	   	
  	 
  	 Pd=den*num; //probability associated to the random number bb
  	 Ptest=rand()/(1.0*RAND_MAX);
  	   	 
  	 if(Ptest<Pd){
  	 	
	   myp[count]=bb;
	   count=count+1;   
  	 	
  	 }
  	 
  	   	
  }
  
  

return myp;
}
//function that fills the initial velocities matrix

void initiate_velocities(matrix& m, double* myp){
	
	//generate random numbers for theta and phi	
	double theta,phi,px,py,pz;	
	int dimr; //dimensions of rows and columns of matrix
	dimr=m.nrows();
	
	for(int i=0;i<dimr;i++){
		//generate random angles
		theta=(rand()/(1.0*RAND_MAX))*pi;
	    phi=(rand()/(1.0*RAND_MAX))*(2*pi); 
	    
	    //variables to fill the matrix
	    px=myp[i]*(sin(theta)*cos(phi));
	    py=myp[i]*(sin(theta)*sin(phi));
	    pz=myp[i]*cos(theta); 
		//components one by one, filled with the appropiate value
		//x component
		m.fill(px,i,0);
		//y component
		m.fill(py,i,1);	  
		//z component
		m.fill(pz,i,2); 
		
	}
	
	//ERASE THIS LINE WHEN PROGRAM IS COMPLETED
	//m.printm();
	
}

//generate position of particles
void initiate_positions(matrix& m,int Npart){
	int count,ylim,zlim;
	double density,epsilon,dx,x,y,z;//density and minimum radius associated to each particle, and epsilon is the value in order to avoid one particle lying on one of
	//the faces of the cube;
	epsilon=0.1;
	//The box lives in the interval -1/2 + 1/2 for each coordinate.
	density=1.0/(Npart);
    //We need to generate NxNxN positions (a cube). We start by creating the points in the z=+epsilon level, and we rise them till z= 1-epsilon. 
    dx=0.8/(49.0); //element separating two points in the initial positions cube
    
    matrix rowpos(50,1);
    x=0.0;
    
    for(int i=0;i<50;i++){
    	rowpos.fill(x,i,0);
    	x=x+dx;
    }
    
    z=0;//Initialize at z=0 plane
    count=0;
    y=0;
    ylim=49;
    zlim=2499;
    
	for(int k=0;k<Npart;k++){
		
		m.fill(rowpos.get(count,0),k,0);
	    m.fill(y,k,1);//fill y  and z values
	    m.fill(z,k,2);
	    
	    if(k==ylim){
	    	ylim=ylim+50;
	    	y=y+dx;	    	
	    	
	    }
		if(k==zlim){
			zlim=zlim+2500;
			y=0.0;
			z=z+dx;
		}
		count=count+1;
		if(count==50){
			count=0;
		}
		
		
	}    
	
}
//store initial data in .txt file

void store1(matrix& m1,matrix& m2,int Npart){
	
    ofstream dataR("pos_data.txt");//store positions in txt file
    ofstream dataV("vel_data.txt");//store velocities
	
	for(int i=0;i<Npart;i++){
		
	dataR<<m1.get(i,0)<<"  "<<m1.get(i,1)<<"  "<<m1.get(i,2);
	dataR<<endl;
	
	dataV<<m2.get(i,0)<<"  "<<m2.get(i,1)<<"  "<<m2.get(i,2);
	dataV<<endl;
	
	}
	
	
}

//function that gets the matrix of relative distances between particles, therefore giving the forces between each of them
double forces_matrix(matrix& r0,matrix& forces,double gamma,double beta,int Npart){
	
	double modulus,F;//modulus and force
	
	for(int i=0;i<Npart;i++){
		for(int j=0;j<Npart;j++){
			if(i==j){
				
				forces.fill(1000000,i,j); //a particle doesnt interact with itself, so the force with this choice will be zero
			}
			
			modulus=sqrt( pow(r0.get(i,0)-r0.get(j,0),2) +  pow(r0.get(i,1)-r0.get(j,1),2) +pow(r0.get(i,2)-r0.get(j,2),2));
			
					
		}
		
	}
}


int main(){
	
	
	srand(time(NULL));//Initialize seed for random numbers
	//Definition of variables/objects    
	
	int Npart,npt,Nmax; //Total number of particles, total number of time frames, maximum number of particles allowed
	double tf,dt,a,gamma,beta;//Final time, differential of time, temp.parameter MB distribution. The parameters gamma and r0 belong to the
	//Buckinham potential, and must be chosen carefully. 
	gamma=0.2;
	beta=0.02;
	
	double* vpoint;//Pointer to velocities distribution
    
	cout<<"Welcome to the C++ program for Molecular Dynamics simulation. Introduce how many particles you want to work with: ";
	cin>>Npart;
	Nmax=125000; //The maximum allowed size of the cube is 50x50x50
	while(Npart>Nmax){
		cout<<"Sorry, number of particles is too big, please reduce it ";
		cin>>Npart;
	}
	
	cout<<"Now choose the value of the parameter a>0 , associated to the temperature: ";
	cin>>a;
	
	while(a<0){
		
		cout<<"Sorry, but the parameter 'a' needs to be positive, try again: ";
		cin>>a;
	}
	
	cout<<endl<<"Include now the total number of time steps you want to take for the simulation: ";
	cin>>npt;
	
	//The initial positions and velocities of the particles. Lets use a Maxwell-Boltzmann distribution for the velocities, whereas the position of particles 
	//should be occupying vertices of a cube in space.
	//Matrices of initial positions and positions at time t, and for velocities., and the matrix of relative distances between particles, relevant for the 
	//forces acting over each individual particle
	
	matrix r0(Npart,3),rt(Npart,3),v0(Npart,3),vt(Npart,3),forces(Npart,Npart);	
	//Generate the array of random velocities
	vpoint=v_rand_MB(Npart,a);
	//Generate matrix of random velocities components
	initiate_velocities(v0,vpoint);
	cout<<"The initial velocities have been allocated in a matrix. "<<endl<<endl;
	initiate_positions(r0,Npart);
	
	store1(r0,v0,Npart);
	
	
	////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	///////////////////////////////////////////////////////////////////////////////////////////////////////////////
	r0.clean();
	v0.clean();
	rt.clean();
	vt.clean();
	forces.clean();
		
	delete[] vpoint;//Free memory block
	
	
	cout<<"End of the program "<<endl;
	
	return 0;
}

[lastchance]

Do you really need to store the forces (line 222)? That matrix is going to be huge. The force on a particle at a particular time and position is only going to be used once: evaluate when needed.

[zaratustra]

Do you really need to store the forces (line 222)? That matrix is going to be huge. The force on a particle at a particular time and position is only going to be used once: evaluate when needed.

That might have caused the crashed, but to be honest, the matrix was going to be "cleaned" at the end by the clean() statement. To be honest, I tried a 100.000 particle size, but the thing is I tried that before and never had problems (because I wanted to see if all the allowed 125.000 maximum particles really filled the cube at t=0). However I can't remember if that was before defining the matrix forces, or afterwards. If it was before, then the addition of that huge matrix (1.25e^2 x 10^10) might have caused the problem.... Is there any way to solve this or have I definitely lost my memory? Also when I try to run a little program now with 2x2 matrices multiplication, it doesn't work, it crashes exactly when multiplication is taking part. Here is the code of that little program:

#include <iostream>
#include <fstream>
#include <cmath>
#include <stdlib.h>
#include "matrix.h"
#include "time.h"

using namespace std;
//OPERATORS OVERLOAD FOR MATRICES
//ADITION
matrix operator+(matrix& M1,matrix& M2){
	int n1,m1,n2,m2;
	double p;//value to fill the matrix
	
	n1=M1.nrows();
	m1=M1.ncolumns();
	n2=M2.nrows();
	m2=M2.ncolumns();
	
	matrix result(n1,m1);
	
	for(int i=0;i<n1;i++){
		for(int j=0;j<m1;j++){
			
			p=M1.get(i,j)+M2.get(i,j);
			result.fill(p,i,j);
			
		}
	}
		
	return result;
		
}
//SUBSTRACTION
matrix operator-(matrix& M1,matrix& M2){
	int n1,m1,n2,m2;
	double p;//value to fill the matrix
	
	n1=M1.nrows();
	m1=M1.ncolumns();
	n2=M2.nrows();
	m2=M2.ncolumns();
	
		
	matrix result(n1,m1);
	
	for(int i=0;i<n1;i++){
		for(int j=0;j<m1;j++){
			
			p=M1.get(i,j)-M2.get(i,j);
			result.fill(p,i,j);
			
		}
	}
		
	return result;
		
}

//MATRIX MULTIPLICATION
matrix operator*(matrix& M1,matrix& M2){
	int n1,m1,n2,m2;
	
	n1=M1.nrows();
	m1=M1.ncolumns();
	n2=M2.nrows();
	m2=M2.ncolumns();
	matrix result(n1,m2);
	double pp;
	
	for(int i=0;i<n1;i++){
		
		for(int j=0;j<m2;j++){
         
		 pp=0.0;
		 for(int k=0;k<m1;k++){
		 	pp=pp+ M1.get(i,k)*M2.get(k,j);
		 }			
		  
		  result.fill(pp,i,j);	
		}
		
	}	
	return result;
	
}

int main(){
	
	srand(time(NULL));	
	matrix A(2,2);
	matrix B(2,2),C;
   
    A.double_rand(0,1);
    
    B.double_rand(0,1);
    A.printm();
    cout<<endl<<endl;
    B.printm();
    
    C=A*B;
     
     C.printm();
   	// clean matrices objects;
	
  
	A.clean();
	B.clean();
	C.clean();
		
	
			
	
	return 0;
}

In this case, it displays both matrices A and B, but when going to display matrix C, it crashes.

Another question is, does the definition of matrix objects inside functions outside the main function involve some memory leak as well?? Thank you for your replies and I hope you can help me with this I am a bit desperate right now since it seems to me I have entirely lost the memory.

[Shadowwolf]

The moment your program exists, the operating system reclaims all memory used by the program. So you shouldn't worry about permanently losing your RAM, this can't happen. You can see this by opening something like the task manager, you should see the memory being reclaimed the moment your application stops.

Could you show the constructor and destructor for your matrix class, this could help solve the problem as this is usually where you allocate/deallocate memory.

[zaratustra]

The moment your program exists, the operating system reclaims all memory used by the program. So you shouldn't worry about permanently losing your RAM, this can't happen. You can see this by opening something like the task manager, you should see the memory being reclaimed the moment your application stops.

But how does that explain that I can't even multiply 2x2 matrices now? I think the first program just screwed up the memory and now there is no space to allocate matrices whenever they are called.

Could you show the constructor and destructor for your matrix class, this could help solve the problem as this is usually where you allocate/deallocate memory.

Yes here is the .cpp file associated with the header "matrix.h"

#include "matrix.h"
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <time.h>

using namespace std;

//allocation of memory

double** matrix::alloc(int rows,int cols){
	x=rows;
	y=cols;	
	
	double** new_matrix=new double*[x];
	for(int i=0;i<x;i++){
		new_matrix[i]=new double[y];
	}
	return new_matrix;
	
}

//Cleaning memory /deallocate
matrix::~matrix(){
	
	for(int i=0;i<x;i++){
		delete[] the_matrix[i];
	}
	delete[] the_matrix;
}
//clean individual matrix 
void matrix::clean(){
	
	for(int i=0;i<x;i++){
		delete[] the_matrix[i];
	}
	delete[] the_matrix;
	
}

//CONSTRUCTORS
//allocate memory
matrix::matrix(){
	
	the_matrix=alloc(0,0);
	
};
//constructor of the matrix as object
matrix::matrix(int rows,int cols){
	
	the_matrix=alloc(rows,cols);
	
}

/// FUNCTIONS
//Print matrix in the command window
void matrix::printm(){
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			cout<<the_matrix[i][j]<<"   ";
		}
		cout<<endl;
		
	}
	cout<<endl<<endl;
	
}

//get matrix elements as double pointer adress

double** matrix::getmatrix(){
	
	return the_matrix;
}
//transpose matrix
void matrix::trans(){
	
	double** transpose=new double*[y];
	for(int i=0;i<y;i++){
		transpose[i]=new double[x];
	}
		
	for(int i=0;i<y;i++){
		for(int j=0;j<x;j++){
			
			transpose[i][j]=the_matrix[j][i];
		}
		
	}
	
	for(int i=0;i<x;i++){
		delete[] the_matrix[i];
	}
	delete[] the_matrix;
	
	the_matrix=transpose;
	
	int new_cols;
	new_cols=x;
	x=y;
	y=new_cols;
			
}

//Number of rows: return
int matrix::nrows(){
	
	return x;
}

//number of columns return:
int matrix::ncolumns(){
	
	return y;
}
//Get element (i,j)
double matrix::get(int row,int column){
	
	return the_matrix[row][column];
	
}

//POinter to element (i,j) of the matrix
void matrix::fill(double value,int row,int column){

the_matrix[row][column]=value;

}

//Populate matrix by elements from command window
void matrix::DOS_populate(){
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			cout<<"Enter elements from left to right :";
			cin>>the_matrix[i][j];
		}
	}
	
}

//flips values of a whole sector of the matrix
void matrix::flip(matrix& m,int row,int column){
	int frow;//final row
	int fcolumn;//final colum;
	
	frow=m.nrows();
	fcolumn=m.ncolumns();
	
	if ((fcolumn>y)||(frow>x)){
	
	cout<<"Error: the matrix you are trying to insert is too big "<<endl;
	exit(EXIT_FAILURE);
	}
	
	
	for(int i=row;i<frow;i++){
		for(int j=column;j<fcolumn;j++){
			
			the_matrix[i][j]=m.get(i,j);
			
		}
	}
	
	
}

//Populate from external array
void matrix::populate(double** pp){
	
	
	
	
}

//Identity matrix
void matrix::eye(){
	
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			the_matrix[i][j]=0;
			if(i==j){
				the_matrix[i][j]=1;
			}
		}
		
	}
	
		
}

//zeros matrix
void matrix::zeros(){
	
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			the_matrix[i][j]=0;
			
		}
		
	}
	

}

//Trace of matrix

double matrix::trace(){
    
    double tr;
    tr=0.0;
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			tr=tr + the_matrix[i][j];
		}
		
	}
	return tr;
}
//Sum of matrices
matrix matrix::add(matrix& m){
	
	if(x!=m.nrows() || y!=m.ncolumns()){
		
		cout<<"Matrix dimensions must agree"<<endl;
		exit(EXIT_FAILURE);
	}
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			the_matrix[i][j]+=m.get(i,j);
			
		}
	}
	
}
//substraction of matrices
matrix matrix::subs(matrix& m){
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			the_matrix[i][j]-=m.get(i,j);
			
		}
	}
	
}

matrix matrix::multiply(matrix& m){
	
	
	double** producto=new double*[x];
	
	for(int i=0;i<x;i++){
		producto[i]=new double[m.ncolumns()];
	}
	
		
	if(y!=m.nrows()){
		cout<<" Matrix multiplication fails"<<endl;
	}
	double suma;
	
	for(int i=0;i<x;i++){
		for(int j=0;j<m.ncolumns();j++){
			
			suma=0.0;
			for(int k=0;k<m.nrows();k++){
				
				suma=suma+the_matrix[i][k]*(m.get(k,j));
				
			}
			
			producto[i][j]=suma;
			
		}
		
	}
	
	for(int i=0;i<x;i++){
		delete[] the_matrix[i];
	}
	delete[] the_matrix;
	
	the_matrix=producto;
	y=m.ncolumns();
	
	
}

//Random matrix generator in the interval min - max (DOUBLE)

void matrix::double_rand(double min,double max){
	
	srand(time(NULL));
	double interval,bb;
	double num;//the random number generated	
	
	interval=abs((max-min));
	
	for(int i=0;i<x;i++){
		for(int j=0;j<y;j++){
			
			bb=(rand())/(1.0*RAND_MAX);
			the_matrix[i][j]=interval*bb + min;
		}
		
	}
			
}

Thanks for the reply Shadowwolf

[zaratustra]

UPDATE: Ok, trying new things now, when I run the first code (the molecular dynamics one) I get no errors at all, at least for 1000 particles, AFTER HAVING ERASED THE FORCES MATRICES, as @lastchance suggested. Im now going to try what happens when I go to bigger particles number. In principle, the problem allows a maximum of 125 000 to be positioned in a cube...

UPDATE 2: The program seems to run perfectly even with 125 000 particles, and I included a cout message in the destructor to be sure that each matrix created is being destroyed at the end of the program. The message is displayed satisfactory.

Does that mean everything seems to be working fine, and that the problem was just the inclusion of that huge matrix of forces? In theory this is a matrix of 0.8 GB size and my computer has 16 GB of RAM, thanks to all of those who replied guys!

[lastchance]

Glad it is working. For what you appear to be doing you can loop round the particles, calculate the force on it, divide by the mass and store the ACCELERATION of each particle - that is Npart x 3 elements, not Npart x Npart. Ultimately, you can use the acceleration to update velocity and then position.

If you have a large number N of particles then you need your memory requirements to be scaling as N, not as NxN.

Good luck.

[Shadowwolf]

The crash in multiplying 2x2 matrices isn't due to a memory problem, that crash was due to a segmentation fault. The problem is that your class doesn't follow the "rule of three". It basically means that if your class does any dynamic allocation, you should define a custom copy constructor, a custom copy assignment operator and a destructor. You have the destructor, but lack the other two.

This is important in your example, as the lack of the other two cause the crash in your 2x2 example. The problem is the line:

C=A*B

To help understand why this happens it might easiest to break down exactly what happens:

1. Call matrix operator*(matrix&, matrix&): to calculate A*B, we'll name the result temp
2. Call matrix& matrix::operator=(const matrix&): this is the copy assignment operator, the compiler automatically generates one if you don't define one yourself, it does this to execute C=temp
3. Destroy temp, as we no longer use the object. It does this by calling the destructor.

Now the problem is in what the compiler generates automatically. By default, it performs a "shallow copy", meaning for your class member the_matrix, it doesn't copy the entire array, just the pointer to it. When the destructor for temp execute, you delete[] this, but since the compiler copied this pointer into C, the_matrix of C now holds a pointer to deleted memory, the moment you print C, you access the values, causing the crash.

The crash you see when you run the 2x2 matrix is thus not a bad_alloc, but a segmentation fault.

You could create a simple copy assignment operator to do a "deep copy" (copy the entire array, instead of just the pointer), it would probably look similar to this:

matrix& matrix::operator=(const matrix& a)
{
    this->clean();
    double** new_matrix = this->alloc(a.nrows(), a.ncolumns());
    for(int i = 0; i < a.nrows(); i++)
    {
        for(int j = 0; j < a.ncolumns(); ++j)
        {
            new_matrix[i][j] = a.get(i,j);
        }
    }
    this->the_matrix = new_matrix;
    
    return *this;
}

You will also need to create a copy constructor, which looks quite similar to this one (it has the signature matrix::matrix(const matrix&)). Note that to make this code compile, you'll have to think about the const-correctness of your functions, so you probably want to add a const after nrows, ncolumns and get.

Note that if you're using C++11, you also need to worry about the move assignment operator and move constructor, the rule of three will then become the rule of five.

You can read more about this on wikipedia:
https://en.wikipedia.org/wiki/Rule_of_three_(C%2B%2B_programming)

EDIT:
I was typing this when you got your original program to work. You might still want to take this into consideration though, because you still have the potential for a crash when you use the = operator.

[zaratustra]

Thank you very much for the detailed explanation Shadowwolf, I could have not imagined what was causing the crash after making the first programm to work. So in the end it has to do with operator overloading, in particular the "="? I will go carefully through your explanation and try to learn more about this, I just realised how easy is to misunderstand things when you are new creating classes, thanks very much!

[zaratustra]

Glad it is working. For what you appear to be doing you can loop round the particles, calculate the force on it, divide by the mass and store the ACCELERATION of each particle - that is Npart x 3 elements, not Npart x Npart. Ultimately, you can use the acceleration to update velocity and then position.

If you have a large number N of particles then you need your memory requirements to be scaling as N, not as NxN.

Good luck.

Thank you lastchance, in fact my idea is to implement Verlet's integration for calculating the position of particles, so in one of its versions, even the velocities would not be necessary... Thank you for the suggestion!