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

This code fixes the order problem of the S vector and the other matrix, but their is still a sign problem. Using this matrix for A:

1.000000 1.000000 0.975000
1.000000 -1.000000 0.975000
-1.000000 -1.000000 -0.925000
-1.000000 1.000000 -1.025000

The modified code produces:

s:
2.793961 2.000000 0.035791

This is correct!

For V:

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

This is NOT correct! 

The correct answer for V is:

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

U is also wrong:  the program outputs:

-0.493230 -0.512652 -0.493875
-0.506363 0.487019 0.506368
0.480733 0.512652 -0.506047
0.518861 -0.487019 0.493554

The correct U is:

-0.4932    0.5127    0.4939    
   -0.5064   -0.4870   -0.5064   
    0.4807   -0.5127    0.5060   
    0.5189    0.4870   -0.4936   

Note last column missing for both solutions  for U.



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

>>> Brian Gough <bjg@network-theory.co.uk> 09/12/01 06:35PM >>>
Patch below fixes the problem, I'll do more testing on that routine
tomorrow.


Index: svdstep.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/svdstep.c,v
retrieving revision 1.3
diff -u -r1.3 svdstep.c
--- svdstep.c	2001/08/29 15:37:21	1.3
+++ svdstep.c	2001/09/12 22:33:01
@@ -126,14 +126,14 @@
 
   if (d0 == 0.0)
     {
-      /* Eliminate off-diagonal element in [0,f1;0,d1] */
+      /* Eliminate off-diagonal element in [0,f1;0,d1] to make [d,0;0,0] */
 
       create_givens (f0, d1, &c, &s);
 
-      /* compute B <= G^T B */
+      /* compute B <= G^T B X,  where X = [0,1;1,0] */
 
-      gsl_vector_set (f, 0, c * f0 - s * d1);
-      gsl_vector_set (d, 1, s * f0 + c * d1);
+      gsl_vector_set (d, 0, c * f0 - s * d1);
+      gsl_vector_set (f, 1, s * f0 + c * d1);
       
       /* Compute U <= U G */
 
@@ -145,6 +145,10 @@
 	  gsl_matrix_set (U, i, 1, s * Uip + c * Uiq);
 	}
 
+      /* Compute V <= V X */
+
+      gsl_matrix_swap_columns (V, 0, 1);
+
       return;
     }
   else if (d1 == 0.0)
@@ -194,17 +198,24 @@
           gsl_matrix_set (V, i, 1, s * Vip + c * Viq);
         }
       
-      /* Eliminate off-diagonal elements, work with the column that
-         has the largest norm */
+      /* Eliminate off-diagonal elements, bring column with largest
+         norm to first column */
       
-      if (fabs(a11) + fabs(a21) > fabs(a12) + fabs(a22))
+      if (hypot(a11, a21) < hypot(a12,a22))
         {
-          create_givens (a11, a21, &c, &s);
+          double t1, t2;
+
+          /* B <= B X */
+
+          t1 = a11; a11 = a12; a12 = t1;
+          t2 = a21; a21 = a22; a22 = t2;
+
+          /* V <= V X */
+
+          gsl_matrix_swap_columns(V, 0, 1);
         } 
-      else
-        {
-          create_givens (a22, -a12, &c, &s);
-        }
+
+      create_givens (a11, a21, &c, &s);
       
       /* compute B <= G^T B */