Index: trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- trunk/psLib/src/math/psMinimizePolyFit.c	(revision 21183)
+++ trunk/psLib/src/math/psMinimizePolyFit.c	(revision 24086)
@@ -15,5 +15,5 @@
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  *
- *  XXX: psMatrixLUSolve() does not return error codes when the results are NANs.
+ *  XXX: psMatrixLUSolution() does not return error codes when the results are NANs.
  *
  *  XXX: For clip-fit functions, what should we do if the mask is NULL?
@@ -403,57 +403,30 @@
     }
 
+    bool status = false;
     if (USE_GAUSS_JORDAN) {
-        // GaussJordan version
-        if (!psMatrixGJSolve(A, B)) {
-            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-            goto escape_GJ;
-        } else {
-            // the first nTerm entries in B correspond directly to the desired
-            // polynomial coefficients.  this is only true for the 1D case
-            for (psS32 k = 0; k < numTerms; k++) {
-                myPoly->coeff[k] = B->data.F64[k];
-                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
-            }
-            // The constant needs to be multiplied by 2, because it's half the a_0.
-            myPoly->coeff[0] *= 2.0;
-            myPoly->coeffErr[0] *= 2.0;
-        }
+        status = psMatrixGJSolve(A, B);
     } else {
-        // LUD version of the fit
-        psImage *ALUD = NULL;
-        psVector* outPerm = NULL;
-        psVector* coeffs = NULL;
-
-        ALUD = psImageAlloc(numTerms, numTerms, PS_TYPE_F64);
-        ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        if (ALUD == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            goto escape_LUD;
-        } else {
-            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-            if (coeffs == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                goto escape_LUD;
-            } else {
-                for (psS32 k = 0; k < numTerms; k++) {
-                    myPoly->coeff[k] = coeffs->data.F64[k];
-                }
-            }
-        }
-        psFree(ALUD);
-        psFree(coeffs);
-        psFree(outPerm);
-    }
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+	goto escape;
+    } 
+
+    // the first nTerm entries in B correspond directly to the desired
+    // polynomial coefficients.  this is only true for the 1D case
+    for (psS32 k = 0; k < numTerms; k++) {
+	myPoly->coeff[k] = B->data.F64[k];
+	myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
+    }
+    // The constant needs to be multiplied by 2, because it's half the a_0.
+    myPoly->coeff[0] *= 2.0;
+    myPoly->coeffErr[0] *= 2.0;
+
     psFree(A);
     psFree(B);
     return true;
 
-escape_LUD:
-    // XXX drop the LUD version!
-    psFree(A);
-    psFree(B);
-    return false;
-
-escape_GJ:
+escape:
     psFree(A);
     psFree(B);
@@ -628,64 +601,32 @@
     }
 
-    // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
+    bool status = false;
     if (USE_GAUSS_JORDAN) {
-        // GaussJordan version
-        if (!psMatrixGJSolve(A, B)) {
-            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-            goto escape_GJ;
-        } else {
-            // the first nTerm entries in B correspond directly to the desired
-            // polynomial coefficients.  this is only true for the 1D case
-            for (psS32 k = 0; k < nTerm; k++) {
-		if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
-		myPoly->coeff[k] = B->data.F64[k];
-		myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
-            }
-        }
+        status = psMatrixGJSolve(A, B);
     } 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);
-        if (ALUD == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            goto escape_LUD;
-        } else {
-            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-            if (coeffs == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                goto escape_LUD;
-            } else {
-                for (psS32 k = 0; k < nTerm; k++) {
-		    if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
-                    myPoly->coeff[k] = coeffs->data.F64[k];
-                    // XXX LUD does not give inverse of A
-                }
-            }
-        }
-        psFree(ALUD);
-        psFree(coeffs);
-        psFree(outPerm);
-    }
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+	goto escape;
+    } 
+
+    // the first nTerm entries in B correspond directly to the desired
+    // polynomial coefficients.  this is only true for the 1D case
+    for (psS32 k = 0; k < nTerm; k++) {
+	if (coeffMask[k] & PS_POLY_MASK_FIT) continue;
+	myPoly->coeff[k] = B->data.F64[k];
+	myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
+    }
+
     psFree(A);
     psFree(B);
-
-    psTrace("psLib.math", 4, "---- %s() End ----\n", __func__);
     return true;
 
-escape_LUD:
-    psFree(A);
-    psFree(B);
-    return false;
-
-escape_GJ:
+escape:
     psFree(A);
     psFree(B);
     return false;
 }
-
 
 /******************************************************************************
@@ -1129,26 +1070,31 @@
     }
 
-    if (!psMatrixGJSolve(A, B)) {
-        psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
-        psFree(A);
-        psFree(B);
-        return false;
-    }
-
-    // select the appropriate solution entries
+    bool status = false;
+    if (USE_GAUSS_JORDAN) {
+        status = psMatrixGJSolve(A, B);
+    } else {
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+	goto escape;
+    } 
+
+    // select the appropriate solution entries (retain the incoming values if masked on the fit)
     for (int i = 0; i < nTerm; i++) {
-        int l = i / nYterm;         // x index
-        int m = i % nYterm;         // y index
-
-	// retain the incoming values if masked on the fit
-	if (coeffMask[l][m] & PS_POLY_MASK_FIT) continue;
-	myPoly->coeff[l][m] = B->data.F64[i];
-	myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]);
+        int ix = i / nYterm;         // x index
+        int iy = i % nYterm;         // y index
+	if (coeffMask[ix][iy] & PS_POLY_MASK_FIT) continue;
+	myPoly->coeff[ix][iy] = B->data.F64[i];
+	myPoly->coeffErr[ix][iy] = sqrt(A->data.F64[i][i]);
     }
     psFree(A);
     psFree(B);
-
-    psTrace("psLib.math", 6, "---- %s() end ----\n", __func__);
     return true;
+
+escape:
+    psFree (A);
+    psFree (B);
+    return false;
 }
 
@@ -1523,70 +1469,29 @@
 
 
-    // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
+    bool status = false;
     if (USE_GAUSS_JORDAN) {
-        // does the solution in place
-        // The matrices were overflowing, so I switched to LUD.
-        if (!psMatrixGJSolve(A, B)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
-            goto escape_GJ;
-        }
-        // select the appropriate solution entries
-        for (int i = 0; i < nTerm; i++) {
-            int ix = i / nYZterm; // x index
-            int iy = (i % nYZterm) / nZterm; // y index
-            int iz = (i % nYZterm) % nZterm; // z index
-	    if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
-            myPoly->coeff[ix][iy][iz] = B->data.F64[i];
-            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
-        }
+        status = psMatrixGJSolve(A, B);
     } 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);
-        if (ALUD == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            goto escape_LUD;
-        } else {
-            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-            if (coeffs == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                goto escape_LUD;
-            } else {
-                // select the appropriate solution entries
-                for (psS32 ix = 0; ix < nXterm; ix++) {
-                    for (psS32 iy = 0; iy < nYterm; iy++) {
-                        for (psS32 iz = 0; iz < nZterm; iz++) {
-                            psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
-			    if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
-                            myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
-                            // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
-                            // XXX this is wrong: A is not replaced with inverse(A), as for GaussJordan
-                        }
-                    }
-                }
-            }
-        }
-
-        psFree(ALUD);
-        psFree(coeffs);
-        psFree(outPerm);
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+	goto escape;
+    } 
+
+    // select the appropriate solution entries
+    for (int i = 0; i < nTerm; i++) {
+	int ix = i / nYZterm; // x index
+	int iy = (i % nYZterm) / nZterm; // y index
+	int iz = (i % nYZterm) % nZterm; // z index
+	if (coeffMask[ix][iy][iz] & PS_POLY_MASK_FIT) continue;
+	myPoly->coeff[ix][iy][iz] = B->data.F64[i];
+	myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[i][i]);
     }
     psFree(A);
     psFree(B);
-
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
     return true;
 
-escape_LUD:
-    psFree(A);
-    psFree(B);
-    return false;
-
-escape_GJ:
-
+escape:
     psFree(A);
     psFree(B);
@@ -1986,72 +1891,30 @@
     }
 
-    // XXX: rel10_ifa used psMatrixGJSolve().  However, this failed tests.  So, I'm using psMatrixLUD().
+    bool status = false;
     if (USE_GAUSS_JORDAN) {
-        // does the solution in place
-        // The GaussJordan version was overflowing, so I'm using LUD.
-        if (!psMatrixGJSolve(A, B)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
-            goto escape_GJ;
-        }
-
-        // select the appropriate solution entries
-        for (int i = 0; i < nTerm; i++) {
-            int ix = i / nYZTterm; // x index
-            int iy = (i % nYZTterm) / nZTterm; // y index
-            int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
-            int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
-	    if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
-            myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
-            myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
-        }
+        status = psMatrixGJSolve(A, B);
     } 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);
-        if (ALUD == NULL) {
-            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
-            goto escape_LUD;
-        } else {
-            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-            if (coeffs == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
-                goto escape_LUD;
-            } else {
-                // select the appropriate solution entries
-                for (psS32 ix = 0; ix < nXterm; ix++) {
-                    for (psS32 iy = 0; iy < nYterm; iy++) {
-                        for (psS32 iz = 0; iz < nZterm; iz++) {
-                            for (psS32 it = 0; it < nTterm; it++) {
-                                psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-				if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
-                                myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
-                                // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
-                                // XXX this is wrong: LUD does not supply inverse(A)
-                            }
-                        }
-                    }
-                }
-            }
-        }
-        psFree(ALUD);
-        psFree(coeffs);
-        psFree(outPerm);
+        status = psMatrixLUSolve(A, B);
+    }
+    if (!status) {
+	psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n");
+	goto escape;
+    } 
+
+    // select the appropriate solution entries
+    for (int i = 0; i < nTerm; i++) {
+	int ix = i / nYZTterm; // x index
+	int iy = (i % nYZTterm) / nZTterm; // y index
+	int iz = ((i % nYZTterm) % nZTterm) / nTterm; // z index
+	int it = ((i % nYZTterm) % nZTterm) % nTterm; // t index
+	if (coeffMask[ix][iy][iz][it] & PS_POLY_MASK_FIT) continue;
+	myPoly->coeff[ix][iy][iz][it] = B->data.F64[i];
+	myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[i][i]);
     }
     psFree(A);
     psFree(B);
-
-    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
     return true;
 
-escape_LUD:
-    psFree(A);
-    psFree(B);
-    return false;
-
-escape_GJ:
+escape:
     psFree(A);
     psFree(B);