Получение coredump segfault при перезапуске метода Lanczos после 3 итераций
Поэтому я работаю над сходимостью алгоритма Ланцоша. Я реализовал это на языке C, сначала вычисляю как матрицу V, ортонормированную для подпространства Крылова, так и трехдиагональную симметричную матрицу T, после чего вычисляю собственные значения и собственные векторы T и векторы Ритца, чтобы определить Вектор инициализации для следующей итерации, если допуск не достигнут.
В цикле while я проверяю перед переходом к следующей итерации, достигнут ли допуск или нет. Программа компилируется, но когда я ее выполняю, я получаю дамп ядра segfault после 3 итераций, я компилирую с помощью -g, gdb говорит мне, что у меня есть дамп ядра в цикле, который вычисляет матрицу Ритца (или векторы Ритца, которые являются собственными векторами матрицы 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;
}
1 ответ
У вас повреждение памяти. Вы перезаписываете значение "Ritz[0]" в одной точке. Вы можете узнать это, напечатав значение 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, Крылов имеет N (4) элементов, но каждый подмассив имеет только m (2) элементы. Вот где ты выходишь за пределы.
Теперь я вижу цикл, который проверяет цикл for lи увеличивая k пока идет. Тем не менее, я думаю, что вам не хватает k = 1 в начале внутри внешнего цикла while (строка 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 явно неправильно. Если посмотреть на структуру цикла еще раз, это та же ошибка, что и в прошлый раз, но только с другой переменной.
настройка j = 0 внутри начала внешнего цикла while (снова строка 236) это тоже исправляется.
Похоже, это все с точки зрения немедленного и катастрофического повреждения памяти, потому что после этого программа работает нормально.
На менее важной ноте у вас есть ошибка, это должно быть sizeof(double):
for( int i = 0 ; i < N ; i++ )
Ritz[i] = calloc( m , sizeof(double*) );