Changeset 24080

Timestamp:

May 6, 2009, 2:14:25 PM (17 years ago)

Author:

eugene

Message:

add test for excessive pivot growth (ill-conditioned matrix)

File:

• trunk/Ohana/src/libohana/src/gaussj.c (modified) (6 diffs)

trunk/Ohana/src/libohana/src/gaussj.c

@@ -1 +1 @@
 # include <ohana.h>
+# define GROWTHTEST 0
+# define MAX_RANGE 1.0e7
 
 // Gauss-Jordan elimination using full pivots based on Press et al's description.  Substantially
 // reworked for Ohana: major modifications to conform to C indexing, use a boolean to track the
 // completed pivot rows and catch the singular matrix early on.  Also, much cleaner control loops
-// than their implementation.  XXX this really needs to check on round-off errors (see version by
-// William Kahan
+// than their implementation.  (largely based on version by William Kahan)
+
+// MAX_RANGE is used to test for ill-conditioned input matrices.  For an ill-conditioned
+// matrix, one or more of the pivots trends towards zero.  Rather than allow this to go to the
+// numerical precision, I am raising an error if |growth| > 1e8
 int dgaussjordan (double **A, double **B, int N, int M) {
 
@@ -19 +24 @@
   memset (pivot, 0, N*sizeof(int));
 
+  double growth = 1.0;
+
   // determine underflow conditions
   // double underFlow = DBL_MIN;
-# if (0)
+# if (GROWTHTEST)
   double roundTest = 4.0;
   roundTest /= 3.0;
   roundTest -= 1.0;
   double epsilon = fabs(((roundTest+roundTest) - 1.0) + roundTest);
-  double growth = 1.0;
 # endif
 
@@ -57 +63 @@
     }
 
+# if (GROWTHTEST)
+    fprintf (stderr, "maxcol: %d\n", maxcol);
+    fprintf (stderr, "full A matrix:\n");
+    for (row = 0; row < N; row++) {
+        for (col = 0; col < N; col++) {
+            fprintf (stderr, "%10.3e ", A[row][col]);
+        }
+        fprintf (stderr, "\n");
+    }
+    fprintf (stderr, "\n");
+# endif
+
     // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
     if (pivot[maxcol]) goto escape;
@@ -76 +94 @@
     for (col = 0; col < N; col++) A[maxcol][col] *= tmpval;
     for (col = 0; col < M; col++) B[maxcol][col] *= tmpval;
-    // XXX measure the pivot growth and trigger on over/under flow
-    // growth *= tmpval;
-    // fprintf (stderr, "column: %d, growth: %e, epsilon: %e\n", maxcol, growth, epsilon);
+
+    // check for ill-conditioned matrix
+    growth *= tmpval;
+
+    // report the pivot growth
+#   if (GROWTHTEST)
+    fprintf (stderr, "column: %d, maxval : %f, growth: %e, epsilon: %e\n", maxcol, tmpval, growth, epsilon);
+    fprintf (stderr, "A diagonal: ");
+    for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
+    fprintf (stderr, "\n");
+# endif
+
+    if (fabs(growth) > MAX_RANGE) goto escape;
 
     /* adjust the elements above the pivot */
@@ -122 +150 @@
   memset (pivot, 0, N*sizeof(int));
 
+  float growth = 1.0;
+
   // determine underflow conditions
   // float underFlow = FLT_MIN;
-# if (0)
+# if (GROWTHTEST)
   float roundTest = 4.0;
   roundTest /= 3.0;
   roundTest -= 1.0;
   float epsilon = fabs(((roundTest+roundTest) - 1.0) + roundTest);
-  float growth = 1.0;
 # endif
 
@@ -179 +208 @@
     for (col = 0; col < N; col++) A[maxcol][col] *= tmpval;
     for (col = 0; col < M; col++) B[maxcol][col] *= tmpval;
-    // growth *= tmpval;
-    // fprintf (stderr, "column: %d, growth: %e, epsilon: %e\n", maxcol, growth, epsilon);
+
+    // check for ill-conditioned matrix
+    growth *= tmpval;
+
+    // report the pivot growth
+#   if (GROWTHTEST)
+    fprintf (stderr, "column: %d, maxval : %f, growth: %e, epsilon: %e\n", maxcol, tmpval, growth, epsilon);
+    fprintf (stderr, "A diagonal: ");
+    for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
+    fprintf (stderr, "\n");
+# endif
+
+    if (fabs(growth) > MAX_RANGE) goto escape;
 
     /* adjust the elements above the pivot */