Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 5822)
+++ /trunk/psLib/src/math/psMinimize.c	(revision 5823)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.149 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-20 23:53:03 $
+ *  @version $Revision: 1.150 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-21 00:27:18 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -3347,34 +3347,64 @@
     }
 
-    // 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++) {
-                for (int iz = 0; iz < 2*nZterm; iz++) {
-                    psFree(Sums[ix][iy][iz]);
+
+    if (0) {
+        // does the solution in place
+        // XXX: The GaussJordan version was overflowing, so I'm using 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++) {
+                    for (int iz = 0; iz < 2*nZterm; iz++) {
+                        psFree(Sums[ix][iy][iz]);
+                    }
+                    psFree(Sums[ix][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++) {
-                for (int it = 0; it < nTterm; it++) {
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-                    myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
-                    myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                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++) {
+                    for (int it = 0; it < nTterm; it++) {
+                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
+                        myPoly->coeffErr[ix][iy][iz][it] = 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++) {
+                    for (int it = 0; it < nTterm; it++) {
+                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                        myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
+                        myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                    }
+                }
+            }
+        }
+
+        psFree(ALUD);
+        psFree(coeffs);
+        psFree(outPerm);
+
     }
 
@@ -3469,5 +3499,5 @@
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
     if (z->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, z);
+        PS_VECTOR_GEN_F64_FROM_F32(z64, z);
     } else {
         z64 = (psVector *) z;
@@ -3480,5 +3510,5 @@
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
     if (t->type.type != PS_TYPE_F64) {
-        PS_VECTOR_GEN_F64_FROM_F32(y64, t);
+        PS_VECTOR_GEN_F64_FROM_F32(t64, t);
     } else {
         t64 = (psVector *) t;
Index: /trunk/psLib/test/math/tst_psPolyFit4DOrd.c
===================================================================
--- /trunk/psLib/test/math/tst_psPolyFit4DOrd.c	(revision 5823)
+++ /trunk/psLib/test/math/tst_psPolyFit4DOrd.c	(revision 5823)
@@ -0,0 +1,279 @@
+/*****************************************************************************
+This routine must ensure that the psVectorFitPolynomial4D works correctly.
+We create a vectors of data points (x and y), and populate them with the
+values from an arbitrary function setData().  We then call
+psVectorFitPolynomial4D() with a regular polynomial data structure.  We then
+evaluate the polynomial with the coefficients generated above and determine
+if they are within an error tolerance of the expected values.
+     *****************************************************************************/
+#include <stdio.h>
+#include "pslib_strict.h"
+#include "psTest.h"
+#include "psMemory.h"
+#include "psVector.h"
+#include "psImage.h"
+#include "psMinimize.h"
+#include <math.h>
+#define NUM_DATA 10
+#define POLY_ORDER 3
+#define A 1.0
+#define B 2.0
+#define C 3.0
+#define D 4.0
+#define E 5.0
+#define F 4.0
+#define H 3.0
+#define J 3.0
+#define K 1.0
+#define L 5.0
+#define M 4.0
+#define N 3.0
+#define O 2.0
+#define P 1.0
+#define Q 5.0
+
+#define ERROR_TOLERANCE 0.10
+#define FERR 1.0
+#define VERBOSE 0
+
+psF32 setDataF32(psF32 x, psF32 y, psF32 z, psF32 t)
+{
+    if (0) {
+        // Linear case, for testing.
+        return(A + (B * x) + (C * y) + (D * z) + (E * t));
+    } else {
+        return(A + (B * x) + (C * y) + (D * z) + (E * t) +
+               (F * x * x) + (H * y * y) + (J * z * z) + (K * t * t) +
+               (L * x * y) + (M * x * z) + (N * x * t) + (O * y * z) + (P * y * t) + (Q * z * t));
+    }
+}
+
+psF64 setDataF64(psF64 x, psF64 y, psF64 z, psF64 t)
+{
+    if (0) {
+        // Linear case, for testing.
+        return(A + (B * x) + (C * y) + (D * z) + (E * t));
+    } else {
+        return(A + (B * x) + (C * y) + (D * z) + (E * t) +
+               (F * x * x) + (H * y * y) + (J * z * z) + (K * t * t) +
+               (L * x * y) + (M * x * z) + (N * x * t) + (O * y * z) + (P * y * t) + (Q * z * t));
+    }
+}
+
+psS32 t00F32()
+{
+    psS32 currentId = psMemGetId();
+    psS32 testStatus = true;
+    psS32 memLeaks = 0;
+    psF64 expectData;
+    psF64 actualData;
+
+    psPolynomial4D *myPoly = psPolynomial4DAlloc(POLY_ORDER, POLY_ORDER, POLY_ORDER, POLY_ORDER, PS_POLYNOMIAL_ORD);
+    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+    psVector *z = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+    psVector *t = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+    psVector *f = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+    psVector *fErr = psVectorAlloc(NUM_DATA, PS_TYPE_F32);
+
+    for (psS32 i=0;i<NUM_DATA;i++) {
+        x->data.F32[i] = (psF32) i;
+        y->data.F32[i] = (psF32) i;
+        z->data.F32[i] = (psF32) i;
+        t->data.F32[i] = (psF32) i;
+        f->data.F32[i] = setDataF32(x->data.F32[i], y->data.F32[i], z->data.F32[i], t->data.F32[i]);
+        fErr->data.F32[i] = FERR;
+    }
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial4D(): psF32 data");
+
+    psVectorFitPolynomial4D(myPoly, NULL, 0, f, fErr, x, y, y, t);
+
+    if (VERBOSE) {
+        for (psS32 i=0;i<POLY_ORDER+1;i++) {
+            for (psS32 j=0;j<POLY_ORDER+1;j++) {
+                for (psS32 k=0;k<POLY_ORDER+1;k++) {
+                    for (psS32 m=0;m<POLY_ORDER+1;m++) {
+                        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i][j][k][m]);
+                    }
+                }
+            }
+        }
+    }
+
+    for (psS32 i=0;i<NUM_DATA;i++) {
+        expectData = setDataF32(x->data.F32[i],  y->data.F32[i], z->data.F32[i], t->data.F32[i]);
+        actualData = psPolynomial4DEval(myPoly, x->data.F32[i], y->data.F32[i], z->data.F32[i], t->data.F32[i]);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("TEST ERROR: Fitted data %d: (%.1f %.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F32[i], y->data.F32[i], actualData, expectData);
+            testStatus = false;
+        }
+        if (VERBOSE) {
+            printf("GOOD: Data point [%d] is %f, should be %f.\n", i, actualData, expectData);
+        }
+    }
+
+    if (0) {
+        printf("Running test with fErr set to NULL.\n");
+        psVectorFitPolynomial4D(myPoly, NULL, 0, f, NULL, x, y, z, t);
+        for (psS32 i=0;i<NUM_DATA;i++) {
+            expectData = setDataF32(x->data.F32[i],  y->data.F32[i], z->data.F32[i], t->data.F32[i]);
+            actualData = psPolynomial4DEval(myPoly, x->data.F32[i], y->data.F32[i], z->data.F32[i], t->data.F32[i]);
+            if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+                printf("TEST ERROR: Fitted data %d: (%.1f %.1f %.1f), expected was (%.1f)\n",
+                       i, x->data.F32[i], y->data.F32[i], actualData, expectData);
+                testStatus = false;
+            }
+            if (VERBOSE) {
+                printf("GOOD: Data point [%d] is %f, should be %f.\n", i, actualData, expectData);
+            }
+        }
+    }
+
+    if (0) {
+        psFree(myPoly);
+        printf("Running test with myPoly set to NULL.  Should generate error.\n");
+        myPoly = psVectorFitPolynomial4D(NULL, NULL, 0, f, fErr, NULL, NULL, NULL, NULL);
+        if (myPoly != NULL) {
+            printf("TEST ERROR: psVectorFitPolynomial4D() returned non-NULL.\n");
+            testStatus = false;
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psFree(z);
+    psFree(t);
+    psFree(f);
+    psFree(fErr);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial4D(): equal differences in variable fErr",
+                testStatus);
+
+    return (!testStatus);
+}
+
+psS32 t00F64()
+{
+    psS32 currentId = psMemGetId();
+    psS32 testStatus = true;
+    psS32 memLeaks = 0;
+    psF64 expectData;
+    psF64 actualData;
+
+    psPolynomial4D *myPoly = psPolynomial4DAlloc(POLY_ORDER, POLY_ORDER, POLY_ORDER, POLY_ORDER, PS_POLYNOMIAL_ORD);
+    psVector *x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    psVector *y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    psVector *z = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    psVector *t = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    psVector *f = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    psVector *fErr = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (psS32 i=0;i<NUM_DATA;i++) {
+        x->data.F64[i] = (psF64) i;
+        y->data.F64[i] = (psF64) i;
+        z->data.F64[i] = (psF64) i;
+        t->data.F64[i] = (psF64) i;
+        f->data.F64[i] = setDataF64(x->data.F64[i], y->data.F64[i], z->data.F64[i], t->data.F64[i]);
+        fErr->data.F64[i] = FERR;
+    }
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial4D(): psF64 data");
+
+    psVectorFitPolynomial4D(myPoly, NULL, 0, f, fErr, x, y, y, t);
+
+    if (VERBOSE) {
+        for (psS32 i=0;i<POLY_ORDER+1;i++) {
+            for (psS32 j=0;j<POLY_ORDER+1;j++) {
+                for (psS32 k=0;k<POLY_ORDER+1;k++) {
+                    for (psS32 m=0;m<POLY_ORDER+1;m++) {
+                        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i][j][k][m]);
+                    }
+                }
+            }
+        }
+    }
+
+    for (psS32 i=0;i<NUM_DATA;i++) {
+        expectData = setDataF64(x->data.F64[i],  y->data.F64[i], z->data.F64[i], t->data.F64[i]);
+        actualData = psPolynomial4DEval(myPoly, x->data.F64[i], y->data.F64[i], z->data.F64[i], t->data.F64[i]);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("TEST ERROR: Fitted data %d: (%.1f %.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], y->data.F64[i], actualData, expectData);
+            testStatus = false;
+        }
+        if (VERBOSE) {
+            printf("GOOD: Data point [%d] is %f, should be %f.\n", i, actualData, expectData);
+        }
+    }
+
+    if (0) {
+        printf("Running test with fErr set to NULL.\n");
+        psVectorFitPolynomial4D(myPoly, NULL, 0, f, NULL, x, y, z, t);
+        for (psS32 i=0;i<NUM_DATA;i++) {
+            expectData = setDataF64(x->data.F64[i],  y->data.F64[i], z->data.F64[i], t->data.F64[i]);
+            actualData = psPolynomial4DEval(myPoly, x->data.F64[i], y->data.F64[i], z->data.F64[i], t->data.F64[i]);
+            if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+                printf("TEST ERROR: Fitted data %d: (%.1f %.1f %.1f), expected was (%.1f)\n",
+                       i, x->data.F64[i], y->data.F64[i], actualData, expectData);
+                testStatus = false;
+            }
+            if (VERBOSE) {
+                printf("GOOD: Data point [%d] is %f, should be %f.\n", i, actualData, expectData);
+            }
+        }
+    }
+
+    if (0) {
+        psFree(myPoly);
+        printf("Running test with myPoly set to NULL.  Should generate error.\n");
+        myPoly = psVectorFitPolynomial4D(NULL, NULL, 0, f, fErr, NULL, NULL, NULL, NULL);
+        if (myPoly != NULL) {
+            printf("TEST ERROR: psVectorFitPolynomial4D() returned non-NULL.\n");
+            testStatus = false;
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psFree(z);
+    psFree(t);
+    psFree(f);
+    psFree(fErr);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial4D(): equal differences in variable fErr",
+                testStatus);
+
+    return (!testStatus);
+}
+
+
+psS32 main()
+{
+    psLogSetFormat("HLNM");
+    t00F32();
+    t00F64();
+}