0

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

在 while 循环中，我在进行下一次迭代之前验证是否达到了容差。该程序编译，但是当我执行它时，我在 3 次迭代后得到一个段错误核心转储，我正在使用 -g 编译，gdb 告诉我我在循环中有一个核心转储，用于计算 Ritz 矩阵（或 Ritz 向量是矩阵 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;
}
4

1 回答 1

0

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

然后要找出它被覆盖的位置，在它被分配后设置一个断点（例如第 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。但是，我认为您在k = 1外部 while 循环（第 236 行）的开头缺少 a 。

修复该问题后，您又遇到了一次崩溃。您在第 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 个元素，因此 for 的值j显然是错误的。再次查看循环的结构，它与上次完全相同的错误，但只是使用了不同的变量。

设置j = 0在外部 while 循环的开头（同样是第 236 行）也修复了这个问题。

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

在不太重要的一点上，您在这里有一个错误，应该是sizeof(double)：

 for( int i = 0 ; i < N ; i++ )  
    Ritz[i] = calloc( m , sizeof(double*) );
于 2016-01-16T03:02:45.757 回答