From mboxrd@z Thu Jan  1 00:00:00 1970
From: "Jim Love" <Jim.Love@asml.com>
To: <gsl-discuss@sources.redhat.com>
Subject: Problem with Singular Value Decomposition Algorithm
Date: Wed, 19 Dec 2001 13:20:00 -0000
Message-id: <sb9f3956.087@wiltonhub.svg.com>
X-SW-Source: 2001/msg00491.html

I have downloaded the latest beta release and the SVD algorithm produces the wrong answers.  It appears that columns are swapped in the output and possibly a sign problem.

Here was my test array:

1 1 0.975
1 -1 0.975
-1 -1 -0.925
-1 1 -1.025

The correct answer for the S vector is:  2.7940 2.0000 0.0358
The gls output was: 2.0000 2.7940 0.0358

The Correct Q matrix is:

-0.7155    0.0256   -0.6981
    0.0183    0.9997    0.0179
   -0.6983   -0.0000    0.7158

The gls output was:

-0.025633 -0.715538 -0.698103
-0.999671 0.018347 0.017900
-0.000000 -0.698332 0.715774

A similar problem was seen in the U matrix.  

Any ideas?  Is this caused by my implementation or is it a real bug?

Here is the source code to the program I made to test the LIB:

#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_matrix.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_blas.h>
#include <C:\Program Files\gsl-0-9-2\gsl\gsl_linalg.h>


int main(void)
{

  gsl_matrix *a = gsl_matrix_alloc (4, 3);
  gsl_matrix *u = gsl_matrix_alloc (4, 4);
  gsl_matrix *v = gsl_matrix_alloc (3, 3);
  gsl_vector *s = gsl_vector_alloc (3);
  gsl_vector *work = gsl_vector_alloc (3);

  gsl_matrix_set_zero(a);

  gsl_matrix_set (a, 0, 0, 1);
  gsl_matrix_set (a, 0, 1, 1);
  gsl_matrix_set (a, 0, 2, 0.975);
  gsl_matrix_set (a, 1, 0, 1);
  gsl_matrix_set (a, 1, 1, -1);
  gsl_matrix_set (a, 1, 2, 0.975);
  gsl_matrix_set (a, 2, 0, -1);
  gsl_matrix_set (a, 2, 1, -1);
  gsl_matrix_set (a, 2, 2, -0.925);
  gsl_matrix_set (a, 3, 0, -1);
  gsl_matrix_set (a, 3, 1, 1);
  gsl_matrix_set (a, 3, 2, -1.025);


  printf ("%f %f %f \n",gsl_matrix_get (a, 0, 0), gsl_matrix_get (a, 0, 1), gsl_matrix_get (a, 0, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 1, 0), gsl_matrix_get (a, 1, 1), gsl_matrix_get (a, 1, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 2, 0), gsl_matrix_get (a, 2, 1), gsl_matrix_get (a, 2, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 3, 0), gsl_matrix_get (a, 3, 1), gsl_matrix_get (a, 3, 2));

  gsl_linalg_SV_decomp (a, v, s, work);

  printf ("\n");
  printf ("%f %f %f \n",gsl_matrix_get (a, 0, 0), gsl_matrix_get (a, 0, 1), gsl_matrix_get (a, 0, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 1, 0), gsl_matrix_get (a, 1, 1), gsl_matrix_get (a, 1, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 2, 0), gsl_matrix_get (a, 2, 1), gsl_matrix_get (a, 2, 2));
  printf ("%f %f %f \n",gsl_matrix_get (a, 3, 0), gsl_matrix_get (a, 3, 1), gsl_matrix_get (a, 3, 2));
  

  printf ("\n");
  printf ("%f %f %f \n",gsl_vector_get (s, 0), gsl_vector_get (s, 1), gsl_vector_get (s, 2));

  printf ("\n");
  printf ("%f %f %f \n",gsl_matrix_get (v, 0, 0), gsl_matrix_get (v, 0, 1), gsl_matrix_get (v, 0, 2));
  printf ("%f %f %f \n",gsl_matrix_get (v, 1, 0), gsl_matrix_get (v, 1, 1), gsl_matrix_get (v, 1, 2));
  printf ("%f %f %f \n",gsl_matrix_get (v, 2, 0), gsl_matrix_get (v, 2, 1), gsl_matrix_get (v, 2, 2));
gsl_matrix_free (a);
  gsl_matrix_free (u);
  gsl_matrix_free (v);
  gsl_vector_free (work);
  gsl_vector_free (s);

	return 0;
}

James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659