Thread: Segmentation Fault in 4 dimensional arrange

Hello,

Im working with a 4 dimensional matrix but I'm having a problem with segmentation fault and I cannot find the problem. Could somebody please give me a hand. This is the code:

Code:

#include <stdio.h>
#include <stdlib.h>
#include<math.h>

/* ----- Data for Grid ----- */
#define xmax 70
#define ymax 50
#define Nmax 30
#define Mmax 30

#define Length 0.8 /*In mts*/
#define Height 0.05 /*In mts*/
#define CFL 0.1
#define k 1.380650*pow(10.0,-23)
#define eo 8.854*pow(10.0,-12)
#define e 1.602*pow(10.0,-19)
#define me 9.10938215*pow(10.0,-31)
#define Pi 3.14159265358979
#define Tei 8000.0



/* --- Velocity Vector for Electrons --- */
void Vel(double* Vx, double* Vy, double Vtermax, double DVx, double DVy)
{
int i,j; 
    for (j=1;j<=2*Nmax;j++){
              Vx[j]=DVx*(j-1)-Vtermax;
              Vx[0]=0.0;
              Vx[2*Nmax+1]=0.0;
                 }
                 
    for (j=1;j<=2*Mmax;j++){
              Vy[j]=DVy*(j-1)-Vtermax;
              Vy[0]=0.0;
              Vy[2*Mmax+1]=0.0;
                 }
  }

/* --- Initial Condition for the Distribution Function --- */
void Init(double**** Fd, double**** gxFd, double****gyFd, double**** gvxFd, double**** gvyFd, double* Vx, double* Vy)
{
   int i, j, p, q;
   
   for(i=0; i<=xmax+1; i++){
        for(j=0; j<=ymax+1; j++){
             for(p=0; p<=2*Nmax+1; p++){
                  for(q=0; q<=2*Mmax+1; q++){
      Fd[i][j][p][q]= (1.0/sqrt(2.0*Pi))*exp(-(me*((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0)); 
      gxFd[i][j][p][q]= 0.0;
      gyFd[i][j][p][q]= 0.0;
      gvxFd[i][j][p][q]= -(Vx[p]*(1.0/sqrt(2.0*Pi))*exp(-((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0));
      gvyFd[i][j][p][q]= -(Vy[q]*(1.0/sqrt(2.0*Pi))*exp(-((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0));

                    }
                  }
                }       
              }
              
   /*Boundary conditions*/
     for(i=0;i<=xmax+1;i++){
        for(j=0;j<=ymax+1;j++){
           for(q=0;q<=2*Mmax+1;q++){ 
   Fd[i][j][0][q]= 0.0;
   
   gxFd[i][j][0][q]= 0.0;
   gyFd[i][j][0][q]= 0.0;
   gvxFd[i][j][0][q]= 0.0;
   gvyFd[i][j][0][q]= 0.0;
   
   Fd[i][j][2*Nmax+1][q]= 0.0;
   
   gxFd[i][j][2*Nmax+1][q]= 0.0;
   gyFd[i][j][2*Nmax+1][q]= 0.0; 
   gvxFd[i][j][2*Nmax+1][q]= 0.0;
   gvyFd[i][j][2*Nmax+1][q]= 0.0;                   
                         }
                     }
                 }
   
   /*Boundary Conditions*/
     for(i=0;i<=xmax+1;i++){
        for(j=0;j<=ymax+1;j++){
           for(p=0;p<=2*Nmax+1;p++){ 
   Fd[i][j][p][0]= 0.0;
   
   gxFd[i][j][p][0]= 0.0;
   gyFd[i][j][p][0]= 0.0;
   gvxFd[i][j][p][0]= 0.0;
   gvyFd[i][j][p][0]= 0.0;
   
   Fd[i][j][p][2*Mmax+1]= 0.0;
   
   gxFd[i][j][p][2*Mmax+1]= 0.0;
   gyFd[i][j][p][2*Mmax+1]= 0.0; 
   gvxFd[i][j][p][2*Mmax+1]= 0.0; 
   gvyFd[i][j][p][2*Mmax+1]= 0.0;                  
                         }
                     }
                 }
   
   /*Boundary Condition*/
         for(i=0;i<=xmax+1;i++){ //Aqui puede haber un error
             for(p=0;p<=2*Nmax+1;p++){
                 for(q=0;q<=2*Mmax+1;q++){
   
   Fd[i][0][p][q]= Fd[i][1][p][2*Mmax+1-q];
   
   gxFd[i][0][p][q]= gxFd[i][1][p][2*Mmax+1-q];
   gyFd[i][0][p][q]= gyFd[i][1][p][2*Mmax+1-q];
   gvxFd[i][0][p][q]= gvxFd[i][1][p][2*Mmax+1-q];
   gvyFd[i][0][p][q]= gvyFd[i][1][p][2*Mmax+1-q];
   
   Fd[i][ymax+1][p][q]= Fd[i][ymax][p][2*Mmax+1-q];
   
   gxFd[i][ymax+1][p][q]= gxFd[i][ymax][p][2*Mmax+1-q];
   gyFd[i][ymax+1][p][q]= gyFd[i][ymax][p][2*Mmax+1-q];
   gvxFd[i][ymax+1][p][q]= gvxFd[i][ymax][p][2*Mmax+1-q];
   gvyFd[i][ymax+1][p][q]= gvyFd[i][ymax][p][2*Mmax+1-q];
   
                         }
                    }
              }
              
     }

int main (void)
{
    int i,j,p,q;
    
    /* ----- Data Velocity ----- */
    double Vx[2*Nmax+2], Vy[2*Mmax+2];
    double Vter=sqrt(2.0*k*Tei/(me));
    double Vtermax=4.0*Vter;
    double DVx=2.0*Vtermax/(2*Nmax-1);
    double DVy=2.0*Vtermax/(2*Mmax-1);
    
    /* --- Malloc definition for Distribution Function and Variables higher than 3D ---*/
    double ****Fd;
    double ****gxFd;
    double ****gyFd;
    double ****gvxFd;
    double ****gvyFd;
    
            Fd = (double****)malloc(sizeof(double***)*ymax+2);
            
            gxFd = (double****)malloc(sizeof(double***)*ymax+2);
            
            gyFd = (double****)malloc(sizeof(double***)*ymax+2);
            
            gvxFd = (double****)malloc(sizeof(double***)*ymax+2);
            
            gvyFd = (double****)malloc(sizeof(double***)*ymax+2);
            
            for(i=0;i<ymax+2;i++){
            Fd[i] = (double***)malloc(sizeof(double**)*xmax+2);
            
            gxFd[i] = (double***)malloc(sizeof(double**)*xmax+2);
            
            gyFd[i] = (double***)malloc(sizeof(double**)*xmax+2);
            
            gvxFd[i] = (double***)malloc(sizeof(double**)*xmax+2);
            
            gvyFd[i] = (double***)malloc(sizeof(double**)*xmax+2);
             
             
            for(j=0;j<xmax+2;j++){
            Fd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gxFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gyFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gvxFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gvyFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            for(p=0;p<2*Nmax+2;p++){
            Fd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gxFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gyFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gvxFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gvyFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            
            for(q=0;q<2*Mmax+2;q++){
            Fd[i][j][p][q] = 0.0;
            
            gxFd[i][j][p][q] = 0.0;
            
            gyFd[i][j][p][q] = 0.0;
            
            gvxFd[i][j][p][q] = 0.0;
            
            gvyFd[i][j][p][q] = 0.0;
            
             }                    
             }
             }
             }
    printf("malloc definition\n");
    Vel(Vx, Vy, Vtermax, DVx, DVy);
    Init(Fd, gxFd, gyFd, gvxFd, gvyFd, Vx, Vy);
    printf("Lleno las matrices\n");
    
    system("PAUSE");
    return 0; 
}

Thank you in advance

Code:

    Fd = (double****)malloc(sizeof(double***)*(ymax+2));
    ...
    for(i=0;i<ymax+2;i++){
        Fd[i] = (double***)malloc(sizeof(double**)*(xmax+2));
        ...

You declare ymax+2 as the size of your primary index (left-most set of [ ]). When you access your arrays elsewhere, you use xmax+2 as the end of your primary dimension. Since they're not the same size, you're ending up out of bounds when you try to access Fd[i][...][...][...].

You also need to put all those +2 bits in parentheses. Multiplication happens before addition, so you are currently ending up with effectively (sizeof(double***)*ymax)+2, which is wrong.

anduril462 Thank you very much. I took both of your comments and it actually worked. I was able to run the program with xmax=100, ymax=100, and also with different sizes.

The thing now is, apparently that is the limit because if I try to run the code with lets say xmax=150 and ymax=150, It gives me segmentation fault. Can that be fixed somehow?

These is the code with the changes I made(basically in the malloc definition)

Code:

#include <stdio.h>
#include <stdlib.h>
#include<math.h>

/* ----- Data for Grid ----- */
#define xmax 100
#define ymax 100
#define Nmax 30
#define Mmax 30

#define Length 0.8 /*In mts*/
#define Height 0.05 /*In mts*/
#define CFL 0.1
#define k 1.380650*pow(10.0,-23)
#define eo 8.854*pow(10.0,-12)
#define e 1.602*pow(10.0,-19)
#define me 9.10938215*pow(10.0,-31)
#define Pi 3.14159265358979
#define Tei 8000.0



/* --- Velocity Vector for Electrons --- */
void Vel(double* Vx, double* Vy, double Vtermax, double DVx, double DVy)
{
int i,j; 
    for (j=1;j<=2*Nmax;j++){
              Vx[j]=DVx*(j-1)-Vtermax;
              Vx[0]=0.0;
              Vx[2*Nmax+1]=0.0;
                 }
                 
    for (j=1;j<=2*Mmax;j++){
              Vy[j]=DVy*(j-1)-Vtermax;
              Vy[0]=0.0;
              Vy[2*Mmax+1]=0.0;
                 }
  }

/* --- Initial Condition for the Distribution Function --- */
void Init(double**** Fd, double**** gxFd, double****gyFd, double**** gvxFd, double**** gvyFd, double* Vx, double* Vy)
{
   int i, j, p, q;
   
   for(i=0; i<=xmax+1; i++){
        for(j=0; j<=ymax+1; j++){
             for(p=0; p<=2*Nmax+1; p++){
                  for(q=0; q<=2*Mmax+1; q++){
      Fd[i][j][p][q]= (1.0/sqrt(2.0*Pi))*exp(-(me*((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0)); 
      gxFd[i][j][p][q]= 0.0;
      gyFd[i][j][p][q]= 0.0;
      gvxFd[i][j][p][q]= -(Vx[p]*(1.0/sqrt(2.0*Pi))*exp(-((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0));
      gvyFd[i][j][p][q]= -(Vy[q]*(1.0/sqrt(2.0*Pi))*exp(-((Vx[p]*Vx[p])+(Vy[q]*Vy[q]))/2.0));

                    }
                  }
                }       
              }
              
   /*Boundary conditions*/
     for(i=0;i<=xmax+1;i++){
        for(j=0;j<=ymax+1;j++){
           for(q=0;q<=2*Mmax+1;q++){ 
   Fd[i][j][0][q]= 0.0;
   
   gxFd[i][j][0][q]= 0.0;
   gyFd[i][j][0][q]= 0.0;
   gvxFd[i][j][0][q]= 0.0;
   gvyFd[i][j][0][q]= 0.0;
   
   Fd[i][j][2*Nmax+1][q]= 0.0;
   
   gxFd[i][j][2*Nmax+1][q]= 0.0;
   gyFd[i][j][2*Nmax+1][q]= 0.0; 
   gvxFd[i][j][2*Nmax+1][q]= 0.0;
   gvyFd[i][j][2*Nmax+1][q]= 0.0;                   
                         }
                     }
                 }
   
   /*Boundary Conditions*/
     for(i=0;i<=xmax+1;i++){
        for(j=0;j<=ymax+1;j++){
           for(p=0;p<=2*Nmax+1;p++){ 
   Fd[i][j][p][0]= 0.0;
   
   gxFd[i][j][p][0]= 0.0;
   gyFd[i][j][p][0]= 0.0;
   gvxFd[i][j][p][0]= 0.0;
   gvyFd[i][j][p][0]= 0.0;
   
   Fd[i][j][p][2*Mmax+1]= 0.0;
   
   gxFd[i][j][p][2*Mmax+1]= 0.0;
   gyFd[i][j][p][2*Mmax+1]= 0.0; 
   gvxFd[i][j][p][2*Mmax+1]= 0.0; 
   gvyFd[i][j][p][2*Mmax+1]= 0.0;                  
                         }
                     }
                 }
   
   /*Boundary Condition*/
         for(i=0;i<=xmax+1;i++){ //Aqui puede haber un error
             for(p=0;p<=2*Nmax+1;p++){
                 for(q=0;q<=2*Mmax+1;q++){
   
   Fd[i][0][p][q]= Fd[i][1][p][2*Mmax+1-q];
   
   gxFd[i][0][p][q]= gxFd[i][1][p][2*Mmax+1-q];
   gyFd[i][0][p][q]= gyFd[i][1][p][2*Mmax+1-q];
   gvxFd[i][0][p][q]= gvxFd[i][1][p][2*Mmax+1-q];
   gvyFd[i][0][p][q]= gvyFd[i][1][p][2*Mmax+1-q];
   
   Fd[i][ymax+1][p][q]= Fd[i][ymax][p][2*Mmax+1-q];
   
   gxFd[i][ymax+1][p][q]= gxFd[i][ymax][p][2*Mmax+1-q];
   gyFd[i][ymax+1][p][q]= gyFd[i][ymax][p][2*Mmax+1-q];
   gvxFd[i][ymax+1][p][q]= gvxFd[i][ymax][p][2*Mmax+1-q];
   gvyFd[i][ymax+1][p][q]= gvyFd[i][ymax][p][2*Mmax+1-q];
   
                         }
                    }
              }
              
     }

int main (void)
{
    int i,j,p,q;
    
    /* ----- Data Velocity ----- */
    double Vx[2*Nmax+2], Vy[2*Mmax+2];
    double Vter=sqrt(2.0*k*Tei/(me));
    double Vtermax=4.0*Vter;
    double DVx=2.0*Vtermax/(2*Nmax-1);
    double DVy=2.0*Vtermax/(2*Mmax-1);
   
   int xmall, ymall, Nmall, Mmall;
   xmall=xmax+2;
   ymall=ymax+2;
   Nmall=(2*Nmax)+2;
   Mmall=(2*Mmax)+2;
    
    /* --- Malloc definition for Distribution Function and Variables higher than 3D ---*/
    double ****Fd;
    double ****gxFd;
    double ****gyFd;
    double ****gvxFd;
    double ****gvyFd;
    
            Fd = (double****)malloc(sizeof(double***)*(xmax+2));
            
            gxFd = (double****)malloc(sizeof(double***)*(xmax+2));
            
            gyFd = (double****)malloc(sizeof(double***)*(xmax+2));
            
            gvxFd = (double****)malloc(sizeof(double***)*(xmax+2));
            
            gvyFd = (double****)malloc(sizeof(double***)*(xmax+2));
            
            for(i=0;i<xmax+2;i++){
            Fd[i] = (double***)malloc(sizeof(double**)*(ymax+2));
            
            gxFd[i] = (double***)malloc(sizeof(double**)*(ymax+2));
            
            gyFd[i] = (double***)malloc(sizeof(double**)*(ymax+2));
            
            gvxFd[i] = (double***)malloc(sizeof(double**)*(ymax+2));
            
            gvyFd[i] = (double***)malloc(sizeof(double**)*(ymax+2));
             
             
            for(j=0;j<ymax+2;j++){
            Fd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gxFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gyFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gvxFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            gvyFd[i][j] = (double**)malloc(sizeof(double*)*(2*Nmax+2));
            
            for(p=0;p<2*Nmax+2;p++){
            Fd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gxFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gyFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gvxFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            gvyFd[i][j][p] = (double*)malloc(sizeof(double)*(2*Mmax+2));
            
            
            for(q=0;q<2*Mmax+2;q++){
            Fd[i][j][p][q] = 0.0;
            
            gxFd[i][j][p][q] = 0.0;
            
            gyFd[i][j][p][q] = 0.0;
            
            gvxFd[i][j][p][q] = 0.0;
            
            gvyFd[i][j][p][q] = 0.0;
            
             }                    
             }
             }
             }
    printf("malloc definition\n");
    Vel(Vx, Vy, Vtermax, DVx, DVy);
    Init(Fd, gxFd, gyFd, gvxFd, gvyFd, Vx, Vy);
    printf("Lleno las matrices\n");
    
    system("PAUSE");
    return 0; 
}

Its working for the current values of xmax, ymax, Nmax and Mmax. Now, I increase the value of one of those, say xmax=150 (or any other) I get segmentation fault.

I will give a read to the FAQ, thank you

- I keep casting the return of malloc because I didn't know it wasn't recommended until "Salem" wrote about it. That's how initially I learned how to use malloc
- The thing with the "+2" has to do with my boundary conditions which are not important now. Its just to sate that my domain(lets say in x) has a size of xmax and an additional two grids for boundary conditions, that's all.
- The loop that starts in 1 is just in one function and I'm fixing the [0] element and the last as well.
- And I'm getting the segmentation fault exactly in line 199 in the last loop of the malloc in "gxFd[i][j][p][q] = 0.0;"

I'm not so skillful, that's why I asked things that might look obvious

Ok
Thank you so much for your help, It has actually been very useful as it helped me to understand more about pointers. Nevertheless, the problem still persist and I want to know a few thinks I've been wondering

1. You recommended not to cast the return of malloc when using C. I was using DevC++ to compile my programs and also Visual Studio and I get the following error when do that (definition without casting): "cannot convert from void* to double****" etc etc. Now, does this has to do with the compiler?

2. I made an ever better simplified version of the function and instead of using complicated sized I used very simple ones: a=100, b=100, c=100, d=100. The programs works well but if I increase a=200 and b=200 I get segmentation fault.
This is the error I get using Visual Studio: "Unhandled exception at 0x0041382c in Perra.exe: 0xC0000005: Access violation writing location 0x00000000." and it points to the line where I have "P[i][j][k][l]=0.0;"

This is the code (using the casting as I wasnt able to run it without it)

Code:

#include <stdio.h>
#include <stdlib.h>

#define A 200
#define B 200
#define C 100
#define D 100

void main()
{
	double ****P;
	int i,j,k,l;

	printf("A=%d\n",A); printf("B=%d\n",B); printf("C=%d\n",C); printf("D=%d\n",D);

	P=(double****)malloc(A*sizeof(double***));
	for(i=0;i<A;i++)
	{
		P[i]=(double***)malloc(B*sizeof(double**));
		for(j=0;j<B;j++)
		{
			P[i][j]=(double**)malloc(C*sizeof(double*));
			for(k=0;k<C;k++)
			{
				P[i][j][k]=(double*)malloc(D*sizeof(double));
				for(l=0;l<D;l++)
				{
					P[i][j][k][l]=0.0;
				}
			}
		}
	}
//(double****)
	printf("Succesful\n");

	system("pause");
}

Hope I'm not too repetitive but this problem is driving me crazy

Thanks for your help