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

[stefan]

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