Hello all,
I think that there is a bug in gsl_linalg_SV_decomp_jacobi. Given the
following matrix:
1 0 0 0
0 1 0 0
0 0 z 0
0 0 0 1
for z = 0, the singular values computed by SV_decomp_jacobi seem to be wrong:
SV_decomp: [ 1, 1, 1, 0]
SV_decomp_mod: [ 1, 1, 1, 0]
SV_decomp_jacobi: [ 1, 1, 0, 0]
^
WRONG?
for z = 0.1 the singular values are correct but not ordered:
SV_decomp: [ 1, 1, 1, 0.1]
SV_decomp_mod: [ 1, 1, 1, 0.1]
SV_decomp_jacobi: [ 1, 1, 0.1, 1]
I have obtained this result with the program included below using gsl 1.6 and
gsl 1.5. My system is Linux Fedora Core 3 (2.6.10-1.760_FC3). The compilation
options were -Wall `gsl-config --cflags` `gsl-config --libs`
Thank you very much for the GSL!
Alberto Ruiz
/* start of code */
#include <stdio.h>
#include "gsl/gsl_linalg.h"
int
main (void)
{
double z;
z = 0;
/*z = 0.1;*/
double a[] = { 1, 0, 0, 0,
0, 1, 0, 0,
0, 0, z, 0,
0, 0, 0, 1 };
double s[] = { 0, 0, 0, 0};
gsl_matrix_view A = gsl_matrix_view_array(a, 4, 4);
gsl_matrix *U = gsl_matrix_alloc(4,4);
gsl_vector_view S = gsl_vector_view_array(s,4);
gsl_matrix *V = gsl_matrix_alloc(4,4);
gsl_vector *workv = gsl_vector_alloc(4);
gsl_matrix *workm = gsl_matrix_alloc(4,4);
gsl_matrix_memcpy(U,&A.matrix);
gsl_linalg_SV_decomp (U, V, &S.vector, workv);
printf ("SV_decomp: [ %g, %g, %g, %g]\n", s[0], s[1], s[2], s[3]);
gsl_matrix_memcpy(U,&A.matrix);
gsl_linalg_SV_decomp_mod (U, workm, V, &S.vector, workv);
printf ("SV_decomp_mod: [ %g, %g, %g, %g]\n", s[0], s[1], s[2], s[3]);
gsl_matrix_memcpy(U,&A.matrix);
gsl_linalg_SV_decomp_jacobi (U, V, &S.vector);
printf ("SV_decomp_jacobi: [ %g, %g, %g, %g]\n", s[0], s[1], s[2], s[3]);
gsl_vector_free(workv);
gsl_matrix_free(workm);
gsl_matrix_free(V);
gsl_matrix_free(U);
return 0;
}
/* end of code */
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
