[Top][All Lists]

[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 14:59:29 +0100
User-agent: Thunderbird (X11/20081224)

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]