bug-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## [Bug-gsl] bug in file eigen/jacobi.c

 From: stefan Subject: [Bug-gsl] bug in file eigen/jacobi.c Date: Wed, 24 Dec 2008 15:08:32 +0100 User-agent: Thunderbird 2.0.0.18 (X11/20081224)

```Hi
I have found a bug in the function:

```
int gsl_eigen_jacobi (gsl_matrix * a, gsl_vector * eval, gsl_matrix * evec, unsigned int max_rot, unsigned int *nrot)
```
```
The Problem is that the routine does not stop if a good solution is found! That means if I set the parameter max_rot for example to 1000 and i call the function then nrot is set to 1000, so nrot is always max_rot. I figured out the the function:
```inline static double norm (gsl_matrix * A)
is probably wrong.

Stefan Kolb
```
```#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>

int main (void)
{
size_t max_it = 10000, done_it = 0;
size_t N = 3;
double data[] = {1, 2, 3,
4, 5, 6,
7, 8, 9,};

gsl_matrix_view m = gsl_matrix_view_array (data, N, N);
gsl_vector *eval = gsl_vector_alloc (N);
gsl_matrix *evec = gsl_matrix_alloc (N, N);

gsl_eigen_jacobi (&m.matrix, eval, evec,max_it, &done_it);
printf("max iterations\t%i\t done iterations %i\n\n",max_it,done_it);

{
int i;

for (i = 0; i < N; i++)
{
double eval_i = gsl_vector_get (eval, i);
gsl_vector_view evec_i  = gsl_matrix_column (evec, i);

printf ("eigenvalue = %g\n", eval_i);
printf ("eigenvector = \n");
gsl_vector_fprintf (stdout,
&evec_i.vector, "%g");
}
}

gsl_vector_free (eval);
gsl_matrix_free (evec);

return 0;
}
```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]