SVD bugs in gsl version .9.1?

[Devin Balkcom <devin at ri dot cmu dot edu>]

Hi,

I'm using gsl to take singular value decompositions.  I have had
problems with all three of the
implementations.   I have written a test program.

- gsl_linalg_SV_decomp() does not return for the test matrix.
- gsl_linalg_SV_decomp_mod() also does not return.
- gsl_linalg_SV_decomp_jacobi() does return, although I have had it fail
in other circumstances,
    with an incorrect result.  (I will send test code for this case
later.)

compile with:

gcc svd_test.c -lgsl -lgslcblas -lm

I'm using Redhat 6.2, but I've also had similar problems under Redhat
7.2 and Debian 2.2.
Hope you can help, thanks for a great piece of software!

-Devin

/************** svd_test.c *****************/

#include <gsl/gsl_linalg.h>
#include <stdio.h>

void vector_print(char *s, gsl_vector *v);
void matrix_print(char *s, gsl_matrix *M);

int main(char argc, char **argv) {

  gsl_matrix *A, *V, *X;
  gsl_vector *s, *work;


  A = gsl_matrix_alloc(3, 3);

  gsl_matrix_set_zero(A);

  gsl_matrix_set(A, 0, 0, -.212);
  gsl_matrix_set(A, 0, 1, .977);
  gsl_matrix_set(A, 0, 2, .155);

  matrix_print("A before SVD", A);

  V = gsl_matrix_alloc(3, 3);
  s = gsl_vector_alloc(3);

  X = gsl_matrix_alloc(3, 3);
  work = gsl_vector_alloc(3);

  /* does not return:
  gsl_linalg_SV_decomp(A, V, s, work);
  */

  /* also does not return:
  gsl_linalg_SV_decomp_mod(A, X, V, s, work);
  */

  /* this seems to work: */
  gsl_linalg_SV_decomp_jacobi(A, V, s);

  vector_print("s", s);
  matrix_print("A after SVD (U)", A);
  matrix_print("V", V);



  return 0;

}

void matrix_print(char *s, gsl_matrix *M) {
  int i, j;

  printf("%s (%d x %d):\n", s, M->size1, M->size2);
  for(i = 0; i < M->size1; i++) {
    for(j = 0; j < M->size2; j++) {
      printf("%7.3f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");

}

void vector_print(char *s, gsl_vector *v) {
  int i;
  printf("%s (%d-vector):  ", s, v->size);
  for(i = 0; i < v->size; i++) printf("%f ", gsl_vector_get(v, i));
  printf("\n\n");
}