I sent the wrong version of the program before, which will not reproduce
the error. I attach the correctly faulty one to this mail.
Benedicte Jourdain tried the program as well and answered the following
to my mail:
without the diagonal at 1.0, I 've the same problem on the IMAC
(infinite loop) but not on a machine with cygwin.
After analysis I think the problem is in function
gsl_linalg_householder_transform called from gsl_linalg_symmtd_decomp
called from gsl_eigen_symmv.
This function (gsl_linalg_householder_transform) called the function
hypot which is either the standard library function or gsl_hypot
depending of your configuration.
On my IMAC this is the standard library function which is called. you
can verify this in the file config.h with #define HAVE_DECL_HYPOT 1.
I'am not an official maintener of GSL, so i can only give advices.
After ./configure edit config.h to change #define HAVE_DECL_HYPOT 0
before make, make check and make install
Katrin
On Sunday 08 November 2009, Brian Gough wrote:
> At Wed, 4 Nov 2009 15:50:22 +0100,
>
> Katrin Wolff wrote:
> > I have a problem computing the eigensystem of one particular matrix
> > that I chanced upon. gsl_eigen_symmv() appears to be caught in an
> > infinite loop. I attach a minimal program and the matrix producing
> > the error.
> >
> > Version and architecture are:
> > Architecture: amd64
> > Source: gsl
> > Version: 1.11+dfsg-1
>
> Thanks for the bug report. I have tried the test program on x86 and
> was not able to reproduce the error. I am busy at a conference this
> week, but I will try it again on a 64 bit machine when I get back.
#include<iostream>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_eigen.h>
using namespace std;
int main(int argc, char *argv[]){
size_t length = 45;
gsl_vector * evalues = gsl_vector_calloc(length);
gsl_matrix * evectors = gsl_matrix_calloc(length, length);
gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(length);
gsl_matrix * m = gsl_matrix_alloc(length, length);
gsl_vector * v = gsl_vector_calloc(length);
FILE *f=fopen("/home/wolff/testmatrix.dat","r");
gsl_matrix_fscanf(f,m);
fclose(f);
for(size_t i = 0; i < length; i++){
for(size_t j = 0; j < length; j++){
cout << gsl_matrix_get(m,i,j) << " ";
//if(i==j){
// gsl_matrix_set(m,i,j,1.);
//}
}
cout << endl;
}
cout << "compute eigensystem" << endl;
gsl_eigen_symmv(m, evalues, evectors, w);
cout << "computed eigensystem" << endl;
gsl_eigen_symmv_sort(evalues, evectors, GSL_EIGEN_SORT_VAL_DESC);
gsl_matrix_get_col(v, evectors, 0);
for(size_t i = 0; i < length; i++){
cout << gsl_vector_get(v,i) << endl;
}
for(size_t i = 0; i < length; i++){
cout << "eigenvalues" << gsl_vector_get(evalues,i) << endl;
}
}_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
