Re: gsl_linalg_SV_decomp goes into infinite loop on reasonable data

[tlb at trevorblackwell dot com]

A much simpler matrix for which gsl_linalg_SV_decomp goes into an
infinite loop is:

  0  0  0
  0  0  1
  1  1  0

This is a resaonable matrix to decompose; Scilab's svd of it is:

 v  =
      ! - 0.7071068    0.  - 0.7071068 !
      ! - 0.7071068    0.    0.7071068 !
      !   0.           1.    0.        !
            
 s  =
      !   1.4142136    0.    0. !
      !   0.           1.    0. !
      !   0.           0.    0. !
            
 u  =
      !   0.    0.  - 1. !
      !   0.    1.    0. !
      ! - 1.    0.    0. !
            
which is easily verified.

This is gsl-1.1.1 on FreeBSD 4.5, using gcc 2.95.3. I get the same
with and without -O.

Below is a patch for test.c to try SVDs on various 3x3 matrices. It
provokes an infinite loop on about the 15th matrix.

I don't understand the Golub-Reinsch algorithm well enough to figure
out what the problem is.



--- linalg/test.c.orig Mon Nov 19 13:39:32 2001
+++ linalg/test.c      Thu Jun 13 11:43:36 2002
@@ -203,6 +203,7 @@
 gsl_matrix * row12;
 
 gsl_matrix * A22;
+gsl_matrix * A33;
 
 gsl_matrix_complex * c7;
 
@@ -1385,6 +1386,44 @@
       }
   }
 
+  {
+    int i,j;
+    int iter;
+    double carry;
+    double lower = 0, upper = 1;
+
+    for (i=0; i<3; i++) {
+      for (j=0; j<3; j++) {
+        gsl_matrix_set (A33, i,j, lower);
+      }
+    }
+    
+    for (iter=0; ; iter++) {
+      f = test_SV_decomp_dim(A33, 64 * GSL_DBL_EPSILON);
+      gsl_test(f, "  SV_decomp A(3x3)(%g %g %g | %g %g %g | %g %g %g)",
+               gsl_matrix_get (A33, 0,0),gsl_matrix_get (A33, 0,1),gsl_matrix_get (A33, 0,2),
+               gsl_matrix_get (A33, 1,0),gsl_matrix_get (A33, 1,1),gsl_matrix_get (A33, 1,2),
+               gsl_matrix_get (A33, 2,0),gsl_matrix_get (A33, 2,1),gsl_matrix_get (A33, 2,2));
+      
+      /* increment */
+      carry=1.0;
+      for (i=2; i>=0; i--) {
+        for (j=2; j>=0; j--) {
+          double v=gsl_matrix_get (A33, i,j)+carry;
+          if (v>upper) {
+            v=lower;
+            carry=1.0;
+          } else {
+            carry=0.0;
+          }
+          gsl_matrix_set (A33, i,j, v);
+        }
+      }
+      if (carry) break;
+    }
+    
+  }
+
   return s;
 }
 
@@ -1790,6 +1829,7 @@
   row12 = create_row_matrix(12,12);
 
   A22 = create_2x2_matrix (0.0, 0.0, 0.0, 0.0);
+  A33 = gsl_matrix_alloc(3,3);
 
   /* Matmult now obsolete */
 #ifdef MATMULT


--
Trevor Blackwell       tlb@trevorblackwell.com        (650) 776-7870