# gsl_linalg_SV_decomp goes into infinite loop on reasonable data

tlb@trevorblackwell.com tlb@trevorblackwell.com
Thu Jun 13 12:43:00 GMT 2002

```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

```