Changeset 5823

Timestamp:

Dec 20, 2005, 2:27:18 PM (21 years ago)

Author:

gusciora

Message:

I've incorporated the 4D polynomial fitting routines from EAM into the CVS
tree and added fairly simple tests. I changed from GaussJordan to LUD since
the former was overflowing.

Location:

trunk/psLib

Files:

• src/math/psMinimize.c (modified) (4 diffs)
• test/math/tst_psPolyFit4DOrd.c (added)

trunk/psLib/src/math/psMinimize.c

@@ -10 +10 @@
  *  @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 +3347 @@
     }
 
-    // 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 +3499 @@
     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 +3510 @@
     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;