segmentation fault.memory error when looping over the function

[Kap]

Dear all,
It is been several days that I am stucked with the segmentation error and I do not know anymore what to do I try to do a loop over my function but it is working only for my first iteration. Here the code:

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>
#include <cstring>
#include <sstream>
#include <string>
#include <vector>
#include <gmm/gmm.h>
#include <time.h>

using namespace std;

int main()

{  
int static dim=26;
double *h, *k, *p, *qx, *qy, *c, *r, *dh, *ht;
int N,nx,ny,nt,Ntime,x;
double dx,dy;
dx=0.04;
Ntime=1;
dy=dx;
nx=dim;
ny=dim;
N  = (nx+1)*ny;					
h  = (double *) calloc(N, sizeof(double));
k  = (double *) calloc(N, sizeof(double));
//p  = (double *) calloc(N, sizeof(double));
c  = (double *) calloc(N, sizeof(double));
r  = (double *) calloc(N, sizeof(double));
qx = (double *) calloc(N, sizeof(double));
qy = (double *) calloc(N, sizeof(double));
dh = (double *) calloc(N, sizeof(double));
ht = (double *) calloc(N, sizeof(double));
nt=0;

void solvep(double *h, double *qx, double *qy, int nt,int dim);	
for(int nt=0;nt<=Ntime;nt++)
{		
solvep(h,qx,qy,nt,dim);

	for (x = 0; x <=N; x++)
		{	
		cout<<h[x]<<','<<qx[x]<<','<<qy[x]<<','<<nt<<endl;
		}

}

return 0;
}
void solvep(double *h, double *qx, double *qy, int nt,int dim)

{
	
			

double H,*p;
double dx,dx2,dy;
int    nn, x, y;
   
int static nx,ny;
double h3, qy0, qavg, q,x1,y1;
nx=dim;
ny=dim;
int N;
N=dim*(dim+1);
dx=0.04;
dy=dx;
dx2=dx*dx;
if (nt==0){
p=(double *) calloc(N, sizeof(double));

}
//------------------DEFINE MATRIX TYPE------------------------------------------------------
//------------------Creation of matrix to write and to read---------------------------------
//------------------------------------------------------------------------------------------
//---------Row matrix to read---------------------------------------------------------------
typedef gmm::csr_matrix<double,1> SPMatrix;
//---------To read vector-------------------------------------------------------------------
typedef gmm::rsvector<double> SPVector;
typedef std::vector<double> DVector;
//-------------------------------------------------------------------------------------------
//To read matrix-----------------------------------------------------------------------------
SPMatrix *NN2; 
SPVector *BB2;
//Solution of the LU SYSTEM--------------------------------------------------------------------
DVector *XX2;
//-----------------------Parameters-------------------------------------------------------------
int xx[4],yy[4];
//int dim=25;//nx,ny being the number of nodes in the x and y direction//
//int N=dim*(dim+1);
int nr,ncn,nnz,m,nnz0,itr,n,np,nc,MAXIT;
//int nx,ny,x,y;

double *c,*r;

double annz;
//ONLY Read matrix and vector----------------------------------------------------------------------
//NOT MODIFIABLE--------------------------------------------------------------------------------
NN2 = new SPMatrix(N,N);
BB2 = new SPVector(N);
XX2 = new DVector(N);
//--------Temporary writeable sparse objects:-----------------------------------------------
gmm::row_matrix< gmm::wsvector<double> > NN(N,N);
gmm::wsvector<double> BB(N);
gmm::wsvector<double> XX(N);
xx[0]=-1;
xx[1]=0;	
xx[2]=0;
xx[3]=1;
yy[0]=0;
yy[1]=-1;	
yy[2]=1;
yy[3]=0;
if (nt == 0){
	for (x = 0; x <= nx; x++){
	for (y = 0; y <  ny; y++){
		//H = HMAX*(1 - x*dx/L )- (y*dx-ny*dx/2)*(y*dx-ny*dx/2)/((W*W));
		x1=x*dx;
		y1=y*dy;
		
		//h[x*ny+y] =H > 0 ? H : 0.01;
		if ((x1<=0.4) && (y1>=0.4) &&(y1<=0.6))
		{
		h[x*ny+y] =1.;
		}	
		//cout<<h[x*ny+y]<<endl;
		else 
		{
			h[x*ny+y]=0.2;
		}
						}
						}
					}



/*  Matrix and RHS setup  */
	nnz = 0;
	for (x = 0; x <= nx; x++){						/* Loop over grid */
	for (y = 0; y <  ny; y++){
		nr = x*ny + y;							/* Indexes */
		//cout<<rn[0]<<','<<endl;
		if (x%nx == 0){
			BB[nr]=1.0;	
			NN(nr,nr)=1.0; nnz++;						/
					}
		else{
			NN(nr,nr)=0.0;	
			nnz0 = nnz;
			nnz++;
			

			for (m = 0; m < 4; m++)
			{
				nc = (x+xx[m])*ny + (y+yy[m]+ny)%ny;
				//if (nc>>N){
				//cout<<x<<','<<y<<','<<nc<<','<<endl;}
				annz = 0.5*(h[nr]*h[nr]*h[nr] + h[nc]*h[nc]*h[nc]);
				annz/= 0.5*(2);
				NN(nr,nr)+=  annz;
				NN(nr,nc)=-annz;
				nnz++;
								
			}
		}

		if (x == nx)	BB[nr] =nx;			//1/* Set RHS vector */
		else{
			BB[nr] = 0;
			}
			}
			
}
/*  Linear system solve  */
//linsol(p, N, a, rn, cn, nnz, itr);		
gmm::copy(NN, *NN2);
//cout<<"NN2"<<*NN2;
gmm::clean(NN, 1E-12);
gmm::copy(BB, *BB2);
gmm::clean(BB, 1E-12);
//cout<<"B2"<<*BB2;
//////gmm::ilut_precond< gmm::row_matrix< gmm::rsvector<double> > > P(*M2, 10, 1e-4);
gmm::identity_matrix P;
gmm::iteration iter(1E-8); // defines an iteration object, with a max residu of 1E-8
gmm::gmres(*NN2, *XX2, *BB2, P, 50, iter); // execute the GMRES algorithm <<"X"<<*X2

for (x = 0; x <=N; x++)
{	
p[x]=XX2[0][x];
//cout <<x<<','<<XX2[0][x]<<endl;
}
gmm::clean(*XX2, 1E-12);

for (x = 0; x < nx; x++)					/* x flux */
	for (y = 0; y < ny; y++){
		n  = x*ny + y;
		np = n + ny;
		h3 = 0.5*(h[n]*h[n]*h[n] + h[np]*h[np]*h[np]);
		qx[n] = h3*(p[np] - p[n]);
	}
	for (y = 0; y < ny; y++){					/* Outflow bc */
		n  = (nx-1)*ny + y;
		np = n + ny;
		qx[np] = qx[n];
	}
		
	for (x = 0; x <= nx; x++)					/* y flux */
	for (y = 0; y <  ny; y++){
		n  = x*ny + y;
		np = x*ny + (y+1)%ny;
		h3 = 0.5*(h[n]*h[n]*h[n] + h[np]*h[np]*h[np]);
		qy[n] = h3*(p[np] - p[n]);
	}

	for (x = nx; x > 0; x--)					/* Nodal values of qx */
	for (y = 0; y <  ny; y++){
		n  = x*ny + y;
		qx[n] = 0.5*(qx[n] + qx[n-ny]);
	}

	for (x = 0; x <= nx; x++){				/* Nodal values of qy */
		n  = x*ny;
		qy0 = 0.5*(qy[n] + qy[n+ny-1]);
		for (y = ny-1; y > 0; y--){
			np = n + y;
			qy[np] = 0.5*(qy[np] + qy[np-1]);
		}
		qy[n] = qy0;
	}


	
}	
  

Any advice is welcome, thanks a lot. I would like to iterate over solvep (since h is changing at each iteration..I usually add other fct but now I try to debug why the loop is not working.).The first iteration is working fine but not the other iterations. I
Thanks a lot.

[ne555]

http://www.cplusplus.com/forum/general/112111/
indent your code
¿what's `gmm'?
NN2 = new SPMatrix(N,N); justify the use of dynamic memory.

[ac517]

line 29:

h  = (double *) calloc(N, sizeof(double));

line 45:

for (x = 0; x <=N; x++)
{	
    cout << h[x] << ',' << qx[x] << ',' << qy[x] << ',' << nt << endl;
}

You seem to have some places where you are using <= in your for loops when I think you need < instead.

[Kap]

Hi,
Thanks for looking into.
Gmm is a library that is very easy to install : to solve linear system LU=X.
http://download.gna.org/getfem/html/homepage/download.html

Concerning line 45 thanks ac517 but I have still this segmentation fault error.

I have found it but do not understand why it is working now..
I have put on line 73 p=calloc for every timesteps so I have commented l73 and l75..
Thanks

[ac517]

OK, double* p is a local variable to the solvep function. So every time you enter solvep it hits line 60 and p becomes an uninitialized pointer to a double.

Let me take just the relevant parts as an example:

void solvep(int nt) {
   double *p;
   if(nt==0) {
      p = calloc(100,sizeof(double));
   }
   p[10]=42;
}

int main() {
     solvep(0); /*OK since we will enter the if statement and allocate memory
                          for p before using it*/

     solvep(1); /* Kaboom, now we never enter the if statement so p is 
                          essentially a random pointer to some memory location, 
                          probably not even related to your program, so when you 
                          use p you get a segmentation fault. */  

}

So that is why taking out the if(nt==0) makes it work passed the first iteration. What was the reason for having that there? Why did you originally want to only allocate on the first iteration? If p needs to retain the values from previous calls then you could make p a static variable.

I think you will still have off by one problems, like line 191 where you

for (x = 0; x <=N; x++)
{	
     p[x]=XX2[0][x];
}

These can be even more insidious because C++ doesn't do bounds checking and it is far more likely that the the very next memory location after the end of the array is associated with your program so it won't throw a segmentation fault and will just silently overwrite random parts of your program data. The canonical idiom for iterating through N items in an array is

double* p = calloc(N,sizeof(double));

//YES
for(int i=0;i<N;i++)
{
      cout << p[i] << endl;
}

//NO
for(int i=0;i<=N;i++)
{
      cout << p[i] << endl;
}