Index: /trunk/Ohana/src/libohana/src/gaussj.c
===================================================================
--- /trunk/Ohana/src/libohana/src/gaussj.c	(revision 16119)
+++ /trunk/Ohana/src/libohana/src/gaussj.c	(revision 16120)
@@ -11,6 +11,4 @@
   int *rowIndex;
   int *pivot;
-  int maxcol, maxrow;
-  double maxval, tmpval;
   
   int diag, col, row;
@@ -21,16 +19,24 @@
   memset (pivot, 0, N*sizeof(int));
 
-  // init these so gcc does not complain
-  maxval = 0.0;
-  maxrow = maxcol = 0;
-
-  // we loop along the matrix diagonal, but we do not operate on the diagonal elements in
-  // order instead, we are looking for the current max element and operating on that
-  // diagonal element.  this is effectively column pivoting.  row pivoting is perfomed
-  // explicitly 
+  // determine underflow conditions
+  // double underFlow = DBL_MIN;
+# if (0)
+  double roundTest = 4.0;
+  roundTest /= 3.0;
+  roundTest -= 1.0;
+  double epsilon = fabs(((roundTest+roundTest) - 1.0) + roundTest);
+  double growth = 1.0;
+# endif
+
+  // Following the algorithm laid out by Press et al., we loop along the matrix diagonal,
+  // but we do not operate on the diagonal elements in order instead, we are looking for
+  // the current max element and operating on that diagonal element.  this is effectively
+  // column pivoting.  row pivoting is perfomed explicitly.  
 
   for (diag = 0; diag < N; diag++) {
 
-    maxval = 0.0;
+    double maxval = 0.0;
+    int maxrow = 0;
+    int maxcol = 0;
 
     // search for the next pivot
@@ -63,10 +69,14 @@
     colIndex[diag] = maxcol;
     if (A[maxcol][maxcol] == 0.0) goto escape;
+    // XXX Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow
 
     /* rescale by pivot reciprocal */
-    tmpval = 1.0 / A[maxcol][maxcol];
+    double tmpval = 1.0 / A[maxcol][maxcol];
     A[maxcol][maxcol] = 1.0;
     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);
 
     /* adjust the elements above the pivot */
@@ -111,4 +121,14 @@
   ALLOCATE (pivot, int, N);
   memset (pivot, 0, N*sizeof(int));
+
+  // determine underflow conditions
+  // float underFlow = FLT_MIN;
+# if (0)
+  float roundTest = 4.0;
+  roundTest /= 3.0;
+  roundTest -= 1.0;
+  float epsilon = fabs(((roundTest+roundTest) - 1.0) + roundTest);
+  float growth = 1.0;
+# endif
 
   // we loop along the matrix diagonal, but we do not operate on the diagonal elements in
@@ -152,4 +172,5 @@
     colIndex[diag] = maxcol;
     if (A[maxcol][maxcol] == 0.0) goto escape;
+    // XXX Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow
 
     /* rescale by pivot reciprocal */
@@ -158,4 +179,6 @@
     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);
 
     /* adjust the elements above the pivot */