c - getting a segfault coredump when restarting the Lanczos method after 3 iterations -
so i'm working on convergence of lanczos algorithm. implemented in c language, first calculate both matrix v orthonormal associated krylov subspace, , triadiagonal symetric matrix t, after compute eigenvalues , eigen vectors of t , ritz vectors define init vector upcoming iteration if tolerance isn't reached.
in while loop, verify before proceeding upcoming iteration if tolerance reached or not. program compiles, when execute segfault core dump after 3 iterations, i'm compiling -g, gdb tells me have core dump in loop computes ritz matrix (or ritz vectors matrix eigen vectors):
for(int = 0 ; < n ; i++) ritz[i][k] = temp[i]; i think wrote program properly, don't know problem. i'd appreciate on guys! p.s : 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 = 0 ; < 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; (i = 0; < 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( = 0; < 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( = 0; < 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( = 0; < 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( = 0; < 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( = 0; < 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 = 0 ; < n ; i++) v[i] = w[i]; //return v; } return v; } void print_matrix( double ** , int m , int n) { int i, j; for( = 0; < 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 ** = (double **) malloc ( n * sizeof(double*) ); for( int = 0 ; < n ; i++ ) a[i] = (double*) malloc( n * sizeof(double) ); double ** t = (double **) malloc( m * sizeof(double*) ); for( int = 0 ; < m ; i++ ) t[i] = (double *) malloc ( m * sizeof(double) ); double ** krylov = (double **) malloc (n * sizeof(double*) ); for( int = 0 ; < n ; i++ ) krylov[i] = (double *) malloc( m * sizeof(double) ); double ** ritz = calloc( n , sizeof(double*) ); for( int = 0 ; < n ; i++ ) ritz[i] = calloc( m , sizeof(double*) ); double ** eigenvect_matrix = calloc ( m , sizeof(double*)); for( int = 0 ; < 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 * = 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 before lanczos algorithm\n\n"); for(int = 0 ; < n ; i++) { for(int j = 0 ; j < n ; j++) { printf("%lf\t",a[i][j]); } printf("\n"); } v[0] = 1.000000; for(int = 0 ; < n ; i++) init[i] = 0.000000; for(int = 1 ; < n ; i++) v[i] = 0.000000; for(int = 0 ; < 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 = 1 ; < 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 = 0 ; < 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 = 0 ; < 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 = 0 ; < 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 = 0 ; < m ; i++ ) { for(int j = 0 ; j < m ; j++) { printf("%lf\t",a[i*m+j]); } printf("\n"); } */ //printf("\n\n\tdelinearized eigen vetors\n\n"); for(int = 0 ; < 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 = 0 ; < 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 = 0 ; < 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 = 0 ; < 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 = 0 ; < 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 = 0 ; < n ; i++) { for(int j = 0 ; j < m ; j++) { printf("%lf\t",ritz[i][j]); } printf("\n"); }*/ return 0; }
you have memory corruption. overwriting value of "ritz[0]" @ 1 point. can find out printing value of ritz, ritz[0] in gdb once crash occurs.
then figure out it's being overwritten, set breakpoint after it's allocated (line 183 example) , run program it. use watch ritz[0] command tell gdb break when writes location. let gdb continue. reaches loop:
for(int = 1 ; < n ; i++) { krylov[i][k] = v[i]; } where i=3 , k=4. krylov has n (4) elements, each sub-array has m (2) elements. that's you're going out of bounds.
now, can see loop around for-loop checking against l, , increasing k while goes. however, think you're missing k = 1 @ beginning inside outer while loop (line 236).
after fixing that, reach crash. crash on line 121 (somme += v[j] * a[i][j];), i=0 , j=0. looks a[0] has been overwritten here, , printing in gdb again confirms this. since a parameter, need find out variable came from, looking @ stack bt find came test.c:346:
temp = prod_mat_vect(krylov, n, m, eigenvect, m); so a there krylov. looks krylov[0] got overwritten, let's same watchpoints did ritz[0].
gdb lands on line 256 (t[l][j+1] = beta;), l=0 , j=4. again, t's sub-arrays allocated 2 elements, value j wrong. looking @ structure of loop again, it's exact same bug last time, different variable.
setting j = 0 inside start of outer while loop (again line 236) fixes 1 too.
that seems in terms of immediate , catastrophic memory corruption, because program runs fine after this.
on less important note, have mistake here, should sizeof(double):
for( int = 0 ; < n ; i++ ) ritz[i] = calloc( m , sizeof(double*) );
Comments
Post a Comment