Changeset 37643

Timestamp:

Nov 19, 2014, 7:57:16 PM (12 years ago)

Author:

eugene

Message:

adding option to squish the bad pivots

File:

• branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c (modified) (5 diffs)

branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c

@@ -1 +1 @@
 # include <ohana.h>
-# define GROWTHTEST 1
-# define MAX_RANGE 1.0e9
+# define GROWTHTEST 0
+# define VERY_VERBOSE 0
+# define MAX_RANGE 1.0e16
+# define MIN_PIVOT 1.0e-7
 
 // Gauss-Jordan elimination using full pivots based on Press et al's description.  Substantially
@@ -68 +70 @@
 # 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");
+# if (VERY_VERBOSE)
+    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
 # endif
 
@@ -101 +105 @@
     /* rescale by pivot reciprocal */
     double tmpval = 1.0 / A[maxcol][maxcol];
+
+    if (fabs(A[maxcol][maxcol]) < MIN_PIVOT) {
+      // we are ill-conditioned.  set this row & col to 0.0, 1.0 on pivot 
+      fprintf (stderr, "WARNING: elliminating degenerate pivot %lf @ A[%d][%d]\n", A[maxcol][maxcol], maxcol, maxcol);
+      for (col = 0; col < N; col++) A[maxcol][col] = 0.0;
+      for (col = 0; col < M; col++) B[maxcol][col] = 0.0;
+      for (row = 0; row < N; row++) A[row][maxcol] = 0.0;
+      A[maxcol][maxcol] = 1.0;
+      continue;
+    }
+
+
     // XXX why is this here (don't I double count this element?) A[maxcol][maxcol] = 1.0;
     A[maxcol][maxcol] = 1.0;
@@ -112 +128 @@
 #   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");
+# if (VERY_VERBOSE)
+    fprintf (stderr, "A diagonal: ");
+    for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
+    fprintf (stderr, "\n");
+# endif
 # endif
 
@@ -132 +150 @@
   }
 
-# if (GROWTHTEST)
+# if (GROWTHTEST && VERY_VERBOSE)
   fprintf (stderr, "final 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");
+  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