1

I wrote the program below ti implement an algorithm with big arrays but when N is over 820 i have a segmentation fault problem (core dump) because of memory. I use 3 big arrays to implement my code. How can i fix this error?

#define N 400

void upper(float A[N][N], float x[N], float b[N])
{
    int i,j;
    float sum;
    x[N-1] = b[N-1]/A[N-1][N-1];

    for(i=N-2; i>=0; i--){
        sum = b[i];
        for(j=i+1; j<N; j++){
            sum = sum - x[j]*A[i][j];
        }
        x[i] = sum /A[i][i];
    }
}

void lower(float A[N][N], float x[N], float b[N])
{
    int i,j;
    float sum;
    x[0] = b[0]/A[0][0];
    for(i=1; i<N; i++){
        sum = b[i];
        for(j=0; j<i; j++){
            sum = sum - x[j]*A[i][j];
        }
        x[i] = sum /A[i][i];
    }
}


void cholesky(float A[N][N],float L[N][N])
{
    int i,j,k;
    float sum;

    for(i=0; i<N; i++){
        for(j=i-4; j<i; j++){ //j-4 because i have 0 and i want to do less computations
            if(j<0)
                continue;
            sum = A[i][j];
            for(k=i-4; k<i; k++){
                if(k>=0)
                    sum = sum - L[i][k]*L[j][k];
            }
            L[i][j] = sum /L[j][j];
        }
        sum = A[i][i];
        for(k=i-4; k<i; k++){
            if(k>=0)
                sum = sum - L[i][k]*L[i][k];
        }
        L[i][i] = sqrt(sum);
    }   
}

void transpose(float L[N][N], float LT[N][N])
{
    int i,j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            LT[i][j] = L[j][i];
        }
    }
}

void table(float A[N][N], float aii, float aone, float athree)
{
    int i,j;

    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            if(i==j){
                A[i][j] = aii;
            }
            else if((i+1 ==j) || (i-1==j)){
                A[i][j] = aone;
            }
            else if(i+3==j || i-3==j){
                A[i][j] = athree;
            }
            else{
                A[i][j] = 0;
            }
        }
    }
}

void vector(float b[N], float b1, float b2, float all)
{
    int i;
    for(i=0; i<N; i++){
        b[i] = all;
    }
    b[0] = b[N-1] = b1;
    b[1] = b[N-2] = b2;
    b[2] = b[N-3] = b2;
}

int main()
{
    float A[N][N];
    float L[N][N],LT[N][N];
    float x[N];
    float y[N];
    float b[N];
    int i,j;

    table(A,12,-5,1);
    vector(b,4,-1,0);

    table(L,0,0,0);
    cholesky(A,L);
    transpose(L,LT);
    lower(L,y,b);
    upper(LT,x,y);

    for(i=0; i<N; i++){
        printf("%f\n",x[i]);
    }
}
Lee Yaan
  • 547
  • 7
  • 26
  • You define three matrices and three vectors on the stack. There is probably a stack overflow when `N` is large. You can try to allocate large arrays dynamically on the heap with `malloc`. – M Oehm Dec 08 '16 at 14:39
  • your segmentation fault is coming from the `main` function, where you're declaring the arrays. These variables are declared on the stack which is limited space (`ulimit -a` at your shell prompt for stack space size if you are using linux/unix; typically 4k or 8k) if your variables exceed that size then you're going to have to dynamically allocate the space or increase your stack size using the `ulimit -s ` and see if that "pushes" the limit of your code to a higher `N`. Otherwise use `malloc` to allocate memory explicitly – Ahmed Masud Dec 08 '16 at 14:42

2 Answers2

1

You allocate large matrices on the stack in main. When N is 840, each matrix requires about 2.8 MB. Your stack cannot accommodate these matrices and it overflows.

You could tamper with your system's stack settings, but it is usually better to allocate large matrices and vectors on the heap:

float (*A)[N] = malloc(N * sizeof(*A));
float (*L)[N] = malloc(N * sizeof(*L));
float (*LT)[N] = malloc(N * sizeof(*LT));
float *x = malloc(N * sizeof(*x));
float *y = malloc(N * sizeof(*y));
float *b = malloc(N * sizeof(*b));

Ensure that all arrays have been allocated properly before using them:

if ((A && L && LT && x && y && b) == 0) {
    fprintf(stderr, "Allocation failed.\n");
    exit(1);
}

Don't forget to clean up the resources after you have used them at the end of main:

free(A);
free(L);
free(LT);
free(x);
free(y);
free(b);
M Oehm
  • 28,726
  • 3
  • 31
  • 42
0

Your segmentation fault is caused by a stack overflow. You can use a heap instead. A simple way to use a heap is by declaring your arrays in the global scope.

Refer to this answer on explicitly using heaps: https://stackoverflow.com/a/571961/2791719

Community
  • 1
  • 1
nerraruzi
  • 108
  • 8