getting a segfault coredump when restarting the Lanczos method after 3 iterations

Question

So I'm working on the convergence of Lanczos algorithm. I implemented it in C language, I first calculate both the matrix V which is the orthonormal associated to the Krylov subspace, and the Triadiagonal symetric matrix T, after that I compute the eigenvalues and the eigen vectors of T and the Ritz vectors to define the init vector for the upcoming iteration if the tolerance isn't reached.

In the while loop, I verify before proceeding to the upcoming iteration if the tolerance is reached or not. The program compiles, but when I execute it I get a segfault core dump after 3 iterations, I'm compiling with -g, gdb tells me that I have a core dump in the loop that computes the Ritz matrix (or the Ritz vectors that are the matrix A eigen vectors):

for(int i = 0 ; i < N ; i++) 
            Ritz[i][k] = temp[i];

I think I wrote my program properly, but I don't know where is the problem. I'd appreciate your help on this guys! P.S : to compile gcc -g -Wall -std=c99 test.c -o test -llapacke -lm

the code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <errno.h>
#include <float.h>
#include <lapacke.h>

double Absolute_value(double numb)
{
    if(numb < 0)
        return (double)((-1)*(numb));
    else
        return numb;
}

double Forbenius_norm(double ** A, int N)
{
    double forbenius_norm = 0.000000;
    for(int i = 0 ; i < N ; i++)
    {
        for(int j = 0 ; j < N ; j++)
        {
            forbenius_norm += abs(A[i][j])*abs(A[i][j]);
        }
    }   
    forbenius_norm = sqrt(forbenius_norm);
    return forbenius_norm;
}

double Norme_vector(double * v, int N)
{
    double norme= 0.;
    double somme = 0.;
    int i;
    for (i = 0; i < N; i++){
        somme += (v[i] * v[i]);
    }

    norme = sqrt(somme);
    return norme;
}

double * Prod_Scal_Vect(double * v, int t, double e)
{
    double * prod =(double *) malloc ( t* sizeof( double));
    int i;
    for( i = 0; i < t ; i++){
        prod[i] = e * v[i];
    }
    return prod;
}

double * Div_vect_scal(double * v, int t, double e)
{
    double * prod =(double *) malloc ( t* sizeof( double));
    int i;
    for( i = 0; i < t ; i++){
        prod[i] = v[i]/e;
    }
    return prod;
}

double Self_Dot_Prod(double *v,int t,double *w,int s)
{
    double res = 0.;
    int i;

    if( t != s )
        printf("Erreur : deux vecteurs non conformes au produit scalaire\n");
    else
    { 
        for( i = 0; i < t; i++)
        {
            res += v[i] * w[i];
        }

        //return res;
    }
    return res;
}

double * vect_Sub(double *v,int t,double *w,int s)
{
    double  * res = (double *) malloc ( s * sizeof(double));
    int i;
    if( t != s )
        printf("Erreur : deux vecteurs non conformes à l'addition de vecteurs\n");
    else
    { 
        for( i = 0; i < t; i++)
        {
            res[i] = v[i] - w[i];
        }

        //return res;
    }
    return res;
}

double * Prod_Mat_Vect( double ** A, int nbl, int nbc, double * v, int t)
{
    //double * res = calloc ( nbl, sizeof( double));
    double * res =(double *) malloc ( nbl* sizeof( double));    
    double somme = 0.;
    int i,j;

    if(t != nbc)
    {   perror(" erreur 4");    
        return NULL;
    }
    else
    {
        for( i = 0; i < nbl; i++)
        {
            somme = 0.;

            for( j = 0; j < nbc; j++)
            {
                somme += v[j] * A[i][j];
            }
            res[i] = somme;
        }
        return res;
    }
}

double * assign_vect(double * v, int n, double * w, int m)
{
    if( n != m)
        printf("Erreur: La taille des deux vecteurs doit être identique pour l'affectaion\n");
    else
    {
        for(int i = 0 ; i < n ; i++)
            v[i] = w[i];
        //return v;
    }
    return v;
}

void print_matrix( double ** A , int m , int n) {
        int i, j;

    for( i = 0; i < m; i++ ) 
    {
                for( j = 0; j < n; j++ ) 
        {
            printf( " %lf", A[i][j] );
        }
                printf( "\n" );
        }
}

int main (int argc , char ** argv)
{
    int N=4;    
    int m = 2;
    //int lda = m;

    double ** A = (double **) malloc ( N * sizeof(double*) );

    for( int i = 0 ; i < N ; i++ )  
        A[i] = (double*) malloc( N * sizeof(double) );

    double ** T = (double **) malloc( m * sizeof(double*) );

    for( int i = 0 ; i < m ; i++ )  
        T[i] = (double *) malloc ( m * sizeof(double) );

    double ** Krylov = (double **) malloc (N * sizeof(double*) );

    for( int i = 0 ; i < N ; i++ )  
        Krylov[i] = (double *) malloc( m * sizeof(double) );



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

    for( int i = 0 ; i < N ; i++ )  
        Ritz[i] = calloc( m , sizeof(double*) );

    double ** Eigenvect_matrix = calloc ( m , sizeof(double*));

    for( int i = 0 ; i < m ; i++ )
        Eigenvect_matrix[i] = calloc ( m , sizeof(double));

    double * v = (double *) malloc(N * sizeof(double));
    double * init = (double *) malloc (N * sizeof(double));
    double * gamma = (double *) malloc(N * sizeof(double)); 

    double * eigenvect = calloc(m,sizeof(double));
    double * tab = calloc(m,sizeof(double));
    double * temp = calloc(N,sizeof(double));
    double * a = calloc(m*m,sizeof(double));
    double * w = calloc(m,sizeof(double));

    double beta = 0.000000;
    double alpha = 0.000000;
    int j=0; 
    int k=1;
    int info = -1;
    int n=m,lda=m;
    int count = 0;
    double test_tolerance = 999;

    A[0][0]= 9; A[0][1]=  1;A[0][2]= -2;A[0][3]=  1;
    A[1][0]= 1; A[1][1]=  8;A[1][2]= -3;A[1][3]= -2;
    A[2][0]= -2;A[2][1]= -3;A[2][2]=  7;A[2][3]= -1;
    A[3][0]= 1; A[3][1]= -2;A[3][2]= -1;A[3][3]=  6;

    printf("\n\n\tOriginal Matrix A before Lanczos Algorithm\n\n");
    for(int i = 0 ; i < N ; i++)
        {
                for(int j = 0 ; j < N ; j++)
                {
                        printf("%lf\t",A[i][j]);
                }
                printf("\n");
        }


    v[0] = 1.000000;    

    for(int i = 0 ; i < N ; i++)
        init[i] = 0.000000;

    for(int i = 1 ; i < N ; i++)
        v[i] = 0.000000;

    for(int i = 0 ; i < N ; i++)
        Krylov[i][0] = v[i];

    while(test_tolerance > 0.00001 || count != 50)
    {
        printf("iteration: %d\n",count+1);
        for(int l = 0 ; l < m-1 ; l++)
        {

            gamma = Prod_Mat_Vect(A, N, N, v, N);
            alpha = Self_Dot_Prod(v, N, gamma, N);
            gamma = vect_Sub(gamma,N,vect_Sub(Prod_Scal_Vect(v, N,alpha),N,Prod_Scal_Vect(init, N,beta),N),N);
            init = assign_vect(init, N, v, N); 
            beta = Norme_vector(gamma, N);
            v = Div_vect_scal(gamma, N, beta);

            for(int i = 1 ; i < N ; i++)
            {
                Krylov[i][k] = v[i];
            }
            k++;

            T[l][l] = alpha;
            T[l+1][j] = beta;
            T[l][j+1] = beta;
            j++;
        }   

        gamma = Prod_Mat_Vect(A, N, N, v, N);
        alpha = Self_Dot_Prod(v, N, gamma, N);
        T[m-1][m-1] = alpha;
        gamma = vect_Sub(gamma,N,vect_Sub(Prod_Scal_Vect(v, N,alpha),N,Prod_Scal_Vect(init, N,beta),N),N);
        init = assign_vect(init, N, v, N); 
        test_tolerance = Norme_vector(gamma, N);
        //printf("tolerance1 %lf\n",test_tolerance);
        test_tolerance = Absolute_value(test_tolerance);
        //printf("tolerance2 %lf\n",test_tolerance);    

        /*printf("\n\n\tTridiagonal Matrix T\n\n"); 
        for(int i = 0 ; i < m ; i++)
        {
            for(int j = 0 ; j < m ; j++)
            {
                printf("%lf\t",T[i][j]);
            }
            printf("\n");
        }*/
        //printf("\n\n\tLinearied matrix\n\n");
        for(int i = 0 ; i < m ; i++ )
        {
            for(int j = 0 ; j < m ; j++)
            {
                a[i*m+j] = T[i][j];
                //printf("a[%d][%d]%lf ",i,j,a[i*m+j]);
            }
            //printf("\n");
        }
        /*printf("\n\n\tKrylov Matrix \n\n");

        for(int i = 0 ; i < N ; i++)
        {
            for(int j = 0 ; j < m ; j++)
            {
                printf("%lf\t",Krylov[i][j]);
            }
            printf("\n");
        }
        */  
        info = LAPACKE_dsyev( LAPACK_ROW_MAJOR, 'V', 'U', n, a, lda, w );

        /*printf("\n\n\tEigenvectors\n\n");
        for(int i = 0 ; i < m ; i++ )
        {
            for(int j = 0 ; j < m ; j++)
            {
                printf("%lf\t",a[i*m+j]);
            }
            printf("\n");
        }
        */
        //printf("\n\n\tDelinearized a for eigen vetors\n\n");
        for(int i = 0 ; i < m ; i++ )
                {
                        for(int j = 0 ; j < m ; j++)
                        {
                                Eigenvect_matrix[i][j] = a[i*m+j];
                                //printf("%lf\t",Eigenvect_matrix[i][j]);
                        }
                        //printf("\n");
                }

        for(int i = 0 ; i < m ; i++)
            tab[i] = w[i];  
        /*printf("\n\tEigenvalues\n\n");

        for(int j = 0 ; j < m ; j++)
            printf("%lf\t",tab[j]);
        printf("\n");
        */
        int index = 0;

        for(int k = 0 ; k < m ; k++)
        {   
            for(int i = 0 ; i < m ; i++)
            {
                eigenvect[i] = a[index];    
                index++;
            }

            test_tolerance *= Absolute_value(eigenvect[m-1]); 

            /*for(int j = 0 ; j < m ; j++)
                printf("%lf\t",a[j]);*/

            temp = Prod_Mat_Vect(Krylov, N, m, eigenvect, m);

            for(int i = 0 ; i < N ; i++)
            {   
                Ritz[i][k] = temp[i]; 
                printf("\ti=%d ,\tk=%d\n",i,k);
            }       
        }
        printf("tolerance : %lf\n",test_tolerance); 

        count++;
        printf("\n");       
        for(int i = 0 ; i < N ; i++)
        {
            init[i] = Ritz[i][0] + Ritz[i][1];  
            init[i] /= Norme_vector(init,N);    
        }

    }

    /*printf("\n\n\tRitz vectors\n\n");
    for(int i = 0 ; i < N ; i++)
    {
        for(int j = 0 ; j < m ; j++)
        {
            printf("%lf\t",Ritz[i][j]);
        }
        printf("\n");
    }*/
    return 0;
}

Lapack is a linear algebra library, so it has to be installed on your machine, you can get it by sudo apt-get install liblapacke — Adil, Jan 16 '16 at 01:16

score 1 · Answer 1 · answered Jan 16 '16 at 03:02

You have memory corruption. You are overwriting the value of "Ritz[0]" at one point. You can find that out by printing the value of Ritz, then Ritz[0] in GDB once the crash occurs.

Then to figure out where it's being overwritten, set a breakpoint after it's allocated (line 183 for example) and run the program up to it. Then use the watch Ritz[0] command to tell GDB to break when something writes to that location. Let gdb continue. It reaches this loop:

for(int i = 1 ; i < N ; i++)
{
    Krylov[i][k] = v[i];
}

Where i=3 and k=4. Krylov has N (4) elements, but each sub-array only has m (2) elements. That's where you're going out of bounds.

Now, I can see the loop around that for-loop is checking against l, and increasing k while it goes. However, I think you're missing a k = 1 at the beginning inside the outer while loop (line 236).

After fixing that, you reach another crash. You crash on line 121 (somme += v[j] * A[i][j];), with i=0 and j=0. Looks like A[0] has been overwritten here, and printing it in GDB again confirms this. Since A is a parameter, we need to find out what variable it actually came from, so looking at the stack with bt we find it came from test.c:346:

temp = Prod_Mat_Vect(Krylov, N, m, eigenvect, m);

So A there was Krylov. Looks like Krylov[0] got overwritten, so let's do the same with watchpoints as we did with Ritz[0].

GDB lands us on line 256 (T[l][j+1] = beta;), with l=0 and j=4. Again, T's sub-arrays were only allocated with 2 elements, so that value for j is clearly wrong. Looking at the structure of the loop again, it's the exact same bug as last time, but just with a different variable.

Setting j = 0 inside the start of the outer while loop (again line 236) fixes this one too.

That seems to be all in terms of immediate and catastrophic memory corruption, because the program runs fine after this.

On a less important note, you have a mistake here, this should be sizeof(double):

 for( int i = 0 ; i < N ; i++ )  
    Ritz[i] = calloc( m , sizeof(double*) );

getting a segfault coredump when restarting the Lanczos method after 3 iterations

1 Answers1