Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 5624)
+++ 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);