Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 5818)
+++ trunk/psLib/src/math/psMinimize.c	(revision 5819)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.148 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-20 22:41:06 $
+ *  @version $Revision: 1.149 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-20 23:53:03 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -486,5 +486,4 @@
     psVector *b)
 {
-
     int *indxc,*indxr,*ipiv;
     int Nx, icol, irow;
@@ -512,5 +511,5 @@
         for (j = 0; j < Nx; j++) {
             if (!isfinite(matrix[i][j])) {
-                psError(PS_ERR_UNKNOWN, false, "Input matrix contains NaNs.\n");
+                psError(PS_ERR_UNKNOWN, false, "Input matrix contains NaNs: matrix[%d][%d] is %.2f\n", i, j, matrix[i][j]);
                 goto fescape;
             }
@@ -1430,20 +1429,31 @@
     }
 
-    ySum = 1.0;
-    for (int j = 0; j < nYsum; j++) {
-        xSum = ySum;
+    xSum = 1.0;
+    for (int i = 0; i < nXsum; i++) {
+        ySum = xSum;
+        for (int j = 0; j < nYsum; j++) {
+            sums->data.F64[i][j] = ySum;
+            ySum *= y;
+        }
+        xSum *= x;
+    }
+
+    if (0) {
+        printf("--------------------- BuildSums2D(%.2f %.2f) ---------------------\n", x, y);
         for (int i = 0; i < nXsum; i++) {
-            sums->data.F64[i][j] = xSum;
-            xSum *= x;
-        }
-        ySum *= y;
-    }
+            for (int j = 0; j < nYsum; j++) {
+                printf("(%.2f) ", sums->data.F64[i][j]);
+            }
+            printf("\n");
+        }
+    }
+
     return (sums);
 }
 
 /******************************************************************************
-BuildSums3D(sums, x, y, z, nXterm, nYterm, nZterm): this routine calculates the powers of
-input parameter "x", "y", and "z" between 0 and input parameter nXterms*2,
-nYterm*2, and nZterm*2.  The result is returned as a psImage sums.
+BuildSums3D(sums, x, y, z, nXterm, nYterm, nZterm): this routine calculates
+the powers of input parameter "x", "y", and "z" between 0 and input parameter
+nXterms*2, nYterm*2, and nZterm*2.  The result is returned as a 3-D array sums.
  *****************************************************************************/
 static psF64 ***BuildSums3D(
@@ -1477,17 +1487,46 @@
     // careful with this function: there is no size checking and realloc for reuse
 
-    zSum = 1.0;
-    for (int k = 0; k < nZsum; k++) {
-        ySum = zSum;
-        for (int j = 0; j < nYsum; j++) {
-            xSum = ySum;
-            for (int i = 0; i < nXsum; i++) {
-                sums[i][j][k] = xSum;
-                xSum *= x;
-            }
-            ySum *= y;
-        }
-        zSum *= z;
-    }
+    if (1) {
+        zSum = 1.0;
+        for (int k = 0; k < nZsum; k++) {
+            ySum = zSum;
+            for (int j = 0; j < nYsum; j++) {
+                xSum = ySum;
+                for (int i = 0; i < nXsum; i++) {
+                    sums[i][j][k] = xSum;
+                    xSum *= x;
+                }
+                ySum *= y;
+            }
+            zSum *= z;
+        }
+    } else {
+        xSum = 1.0;
+        for (int i = 0; i < nXsum; i++) {
+            ySum = xSum;
+            for (int j = 0; j < nYsum; j++) {
+                zSum = ySum;
+                for (int k = 0; k < nZsum; k++) {
+                    sums[i][j][k] = zSum;
+                    zSum *= z;
+                }
+                ySum *= y;
+            }
+            xSum *= x;
+        }
+    }
+
+    if (0) {
+        printf("--------------------- BuildSums3D(%.2f %.2f %.2f) ---------------------\n", x, y, z);
+        for (int k = 0; k < nXsum; k++) {
+            for (int j = 0; j < nYsum; j++) {
+                for (int i = 0; i < nZsum; i++) {
+                    printf("(%.2f) ", sums[k][j][i]);
+                }
+                printf("\n");
+            }
+        }
+    }
+
     return (sums);
 }
@@ -2760,8 +2799,10 @@
             for (int iy = 0; iy < nYterm; iy++) {
                 for (int iz = 0; iz < nZterm; iz++) {
-                    if (myPoly->mask[ix][iy][iz])
+                    if (myPoly->mask[ix][iy][iz]) {
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
+                    } else {
+                        int nx = ix + iy*nXterm + iz*nXterm*nYterm;
+                        B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
+                    }
                 }
             }
@@ -2808,33 +2849,56 @@
     }
 
-    //    PS_IMAGE_PRINT_F64(A);
-    //    PS_VECTOR_PRINT_F64(B);
-    // does the solution in place
-    if (false == psGaussJordan (A, B)) {
-        psFree(A);
-        psFree(B);
-
-        for (int ix = 0; ix < 2*nXterm; ix++) {
-            for (int iy = 0; iy < 2*nYterm; iy++) {
-                psFree(Sums[ix][iy]);
-            }
-            psFree(Sums[ix]);
-        }
-        psFree(Sums);
-        psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
-        return(NULL);
-    }
-
-    // select the appropriate solution entries
-    for (int ix = 0; ix < nXterm; ix++) {
-        for (int iy = 0; iy < nYterm; iy++) {
-            for (int iz = 0; iz < nZterm; iz++) {
-                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
-                myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
-            }
-        }
-    }
-
+    if (0) {
+        // does the solution in place
+        // The matrices were overflowing, so I switched to LUD.
+        if (false == psGaussJordan (A, B)) {
+            psFree(A);
+            psFree(B);
+
+            for (int ix = 0; ix < 2*nXterm; ix++) {
+                for (int iy = 0; iy < 2*nYterm; iy++) {
+                    psFree(Sums[ix][iy]);
+                }
+                psFree(Sums[ix]);
+            }
+            psFree(Sums);
+            psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n");
+            return(NULL);
+        }
+        // select the appropriate solution entries
+        for (int ix = 0; ix < nXterm; ix++) {
+            for (int iy = 0; iy < nYterm; iy++) {
+                for (int iz = 0; iz < nZterm; iz++) {
+                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                    myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
+                    myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                }
+            }
+        }
+    } 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);
+
+        // select the appropriate solution entries
+        for (int ix = 0; ix < nXterm; ix++) {
+            for (int iy = 0; iy < nYterm; iy++) {
+                for (int iz = 0; iz < nZterm; iz++) {
+                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                    myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
+                    myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                }
+            }
+        }
+
+        psFree(ALUD);
+        psFree(coeffs);
+        psFree(outPerm);
+    }
     psFree(A);
     psFree(B);