c - 3 次迭代后重新启动 Lanczos 方法时出现段错误 coredump

Question

所以我正在研究 Lanczos 算法的收敛。我用 C 语言实现它，首先计算与 Krylov 子空间相关的正交矩阵 V 和三对角对称矩阵 T，然后计算 T 的特征值和特征向量以及 Ritz vector 来定义如果未达到容差，则为即将到来的迭代初始化 vector 。

在 while 循环中，我在进行下一次迭代之前验证是否达到容差。该程序编译，但当我执行它时，我在 3 次迭代后得到一个段错误核心转储，我使用 -g 进行编译，gdb 告诉我，我在计算 Ritz 矩阵(或 Ritz vector )的循环中有一个核心转储是矩阵 A 特征向量):

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

我认为我的程序编写正确，但我不知道问题出在哪里。我很感激你对这些人的帮助!PS:编译 gcc -g -Wall -std=c99 test.c -o test -llapacke -lm

代码:

#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;
}   

Answer 1

你有内存损坏。您在某一时刻覆盖了“Ritz[0]”的值。您可以通过在崩溃发生后在 GDB 中打印 Ritz 和 Ritz[0] 的值来找出这一点。

然后要找出它被覆盖的位置，请在分配后设置一个断点(例如第 183 行)并运行程序直到该断点。然后使用watch Ritz[0]命令告诉GDB在写入该位置时中断。让gdb继续。它到达这个循环:

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

其中i=3和k=4。 Krylov 有 N (4) 个元素，但每个子数组只有 m (2) 个元素。这就是你要出界的地方。

现在，我可以看到围绕该 for 循环的循环正在检查 l，并在运行时增加 k。但是，我认为您在外部 while 循环内的开头缺少一个 k = 1 (第 236 行)。

修复该问题后，您会遇到另一次崩溃。您在第 121 行 (somme += v[j] * A[i][j];) 崩溃，其中 i=0 和 j=0。看起来A[0]在这里被覆盖了，再次在GDB中打印它证实了这一点。由于 A 是一个参数，我们需要找出它实际来自哪个变量，因此用 bt 查看堆栈，我们发现它来自 test.c :346 :

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

所以A有Krylov。看起来 Krylov[0] 被覆盖了，所以让我们对观察点执行与 Ritz[0] 相同的操作。

GDB 将我们带到第 256 行 (T[l][j+1] = beta;)，其中 l=0 和 j=4。同样，T 的子数组仅分配了 2 个元素，因此 j 的值显然是错误的。再次查看循环的结构，这与上次的错误完全相同，只是变量不同。

在外部 while 循环的开头设置 j = 0(同样是第 236 行)也可以修复此问题。

这似乎就是立即和灾难性的内存损坏，因为程序在此之后运行良好。

<小时/>

在一个不太重要的说明上，你在这里有一个错误，这应该是sizeof(double):

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