Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 5654)
+++ /trunk/psLib/src/math/psMinimize.c	(revision 5655)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.146 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-30 02:00:09 $
+ *  @version $Revision: 1.147 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-01 23:36:23 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1756,12 +1756,31 @@
     }
 
-    // does the solution in place
-    psGaussJordan (A, B);
-
-    // the first nTerm entries in B correspond directly to the desired
-    // polynomial coefficients.  this is only true for the 1D case
-    for (int k = 0; k < nTerm; k++) {
-        myPoly->coeff[k] = B->data.F64[k];
-    }
+    // GaussJordan version
+    if (0) {
+        // does the solution in place
+        psGaussJordan (A, B);
+
+        // the first nTerm entries in B correspond directly to the desired
+        // polynomial coefficients.  this is only true for the 1D case
+        for (int k = 0; k < nTerm; k++) {
+            myPoly->coeff[k] = B->data.F64[k];
+        }
+    } else {
+        // LUD version of the fit
+        psImage *ALUD = NULL;
+        psVector* outPerm = NULL;
+        psVector* coeffs = NULL;
+
+        ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+        ALUD = psMatrixLUD(ALUD, &outPerm, A);
+        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+        for (int k = 0; k < nTerm; k++) {
+            myPoly->coeff[k] = coeffs->data.F64[k];
+        }
+        psFree(ALUD);
+        psFree(coeffs);
+        psFree(outPerm);
+    }
+
 
     psFree(A);
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 5654)
+++ /trunk/psLib/src/math/psStats.c	(revision 5655)
@@ -17,6 +17,6 @@
  *
  *
- *  @version $Revision: 1.155 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-30 02:00:09 $
+ *  @version $Revision: 1.156 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-01 23:35:33 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1379,4 +1379,5 @@
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
+    psTrace(__func__, 4, "range (%f, %f) get (%f)\n", rangeLow, rangeHigh, getThisValue);
     PS_ASSERT_POLY_NON_NULL(myPoly, NAN);
     PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN);
@@ -1386,4 +1387,5 @@
     psF32 fLo = psPolynomial1DEval(myPoly, rangeLow);
     psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh);
+    psTrace(__func__, 4, "function at endpoints are (%f, %f) get (%f)\n", fLo, fHi, getThisValue);
     if (!((fLo <= getThisValue) && (fHi >= getThisValue))) {
         psError(PS_ERR_UNKNOWN, true, "The requested y value (%f) does not fall within the specified range (%f to %f)\n", getThisValue, fLo, fHi);
@@ -1473,5 +1475,8 @@
         y->data.F64[1] = yVec->data.F32[binNum];
         y->data.F64[2] = yVec->data.F32[binNum + 1];
-        psTrace(__func__, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        psTrace(__func__, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
+                xVec->data.F32[binNum],
+                xVec->data.F32[binNum+1],
+                xVec->data.F32[binNum+2]);
         psTrace(__func__, 6, "x vec is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
         psTrace(__func__, 6, "y vec is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
@@ -1538,5 +1543,5 @@
         psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
         psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
-        psTrace(__func__, 6, "Fitted y vec is (%.2f %.2f %.2f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));