This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl_linalg_SV_decomp goes into infinite loop on reasonable data
- From: tlb at trevorblackwell dot com
- To: gsl-discuss at sources dot redhat dot com
- Date: Thu, 13 Jun 2002 12:21:39 -0700
- Subject: Re: gsl_linalg_SV_decomp goes into infinite loop on reasonable data
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