Index: trunk/Ohana/src/libohana/src/gaussj.c
===================================================================
--- trunk/Ohana/src/libohana/src/gaussj.c	(revision 24033)
+++ trunk/Ohana/src/libohana/src/gaussj.c	(revision 24080)
@@ -1,9 +1,14 @@
 # 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,12 +24,13 @@
   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,4 +63,16 @@
     }
 
+# 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,7 +94,17 @@
     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,12 +150,13 @@
   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,6 +208,17 @@
     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 */