Re: Problem with Singular Value Decomposition Algorithm

[Brian Gough <bjg@network-theory.co.uk>]

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