[Bug-gsl] gsl_linalg_SV_decomp_jacobi

[Alberto Ruiz]

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 */