obtenir un coredump segfault lors du redémarrage de la méthode Lanczos après 3 itérations

Question

Donc je travaille sur la convergence de l'algorithme de Lanczos. Je l'ai implémenté en langage C, je calcule d'abord la matrice V qui est l'orthonormal associé au sous-espace de Krylov, et la matrice symétrique triadiale T, après que je calcule les valeurs propres et les vecteurs propres de T et les vecteurs de Ritz pour définir le init vecteur pour l'itération à venir si la tolérance n'est pas atteinte.

Dans la boucle while, je vérifie avant de passer à l'itération à venir si la tolérance est atteinte ou non. Le programme compile, mais quand je l'exécute, je reçois un vidage de base de segfault après 3 itérations, je compile avec -g, gdb me dit que j'ai un core dump dans la boucle qui calcule la matrice de Ritz (ou les vecteurs Ritz sont les vecteurs propres de la matrice A):

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

Je pense que j'ai écrit mon programme correctement, mais je ne sais pas où est le problème. J'apprécierais votre aide sur ce gars! P.S: pour compiler gcc -g -Wall std = c99 test.c -o Test -llapacke -lm

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

Answer 1

Vous avez corruption de mémoire. Vous remplacez la valeur de "Ritz [0]" à un moment donné. Vous pouvez trouver cela en imprimant la valeur de Ritz, puis Ritz[0] dans GDB une fois que le crash se produit. Puis, pour savoir où il est écrasé, définissez un point d'arrêt après l'avoir alloué (ligne 183 par exemple) et exécutez le programme. Utilisez ensuite la commande watch Ritz[0] pour indiquer à GDB de se casser lorsque quelque chose écrit à cet emplacement. Laissez gdb continuer. Il atteint cette boucle:

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

Où i=3 et k=4. Krylov a N (4) éléments, mais chaque sous-tableau a seulement m (2) éléments. C'est là que vous allez hors des limites.

Maintenant, je peux voir la boucle autour de cette boucle fortifie contre l, et en augmentant k pendant qu'il va. Cependant, je pense qu'il vous manque un k = 1 au début dans la boucle while externe (ligne 236).

Après avoir corrigé ce, vous atteignez un autre plantage. Vous vous écrasez sur la ligne 121 (somme += v[j] * A[i][j];), avec i=0 et j=0. On dirait que A[0] a été écrasé ici, et l'imprimer à nouveau dans GDB le confirme. Depuis A est un paramètre, il faut savoir quelle variable il est venu en fait de, regarder si à la pile avec bt nous trouvons est venu de test.c:346:

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

Alors A il y avait Krylov. On dirait que Krylov[0] a été écrasé, alors faisons la même chose avec les points de surveillance que nous avons fait avec Ritz[0]. GDB nous envoie sur la ligne 256 (T[l][j+1] = beta;), avec l=0 et j=4. Encore une fois, les sous-tableaux de T ont seulement été alloués avec 2 éléments, de sorte que la valeur pour j est clairement erronée. En regardant à nouveau la structure de la boucle, c'est exactement le même bug que la dernière fois, mais juste avec une variable différente.

Le réglage j = 0 à l'intérieur du début de la boucle while externe (à nouveau la ligne 236) corrige également celui-ci.

Cela semble être tout en termes de corruption de mémoire immédiate et catastrophique, car le programme fonctionne bien après cela.

---

Sur une note moins importante, vous avez une erreur ici, cela devrait être sizeof(double):

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