Changeset 5827 for branches/eam_rel9_b1/psLib/src/math/psMinimize.c

Timestamp:

Dec 22, 2005, 4:09:39 PM (21 years ago)

Author:

magnier

Message:

added optimization to 3D polynomial fit

File:

• branches/eam_rel9_b1/psLib/src/math/psMinimize.c (modified) (11 diffs)

branches/eam_rel9_b1/psLib/src/math/psMinimize.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.147.4.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-19 19:44:27 $
+ *  @version $Revision: 1.147.4.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-23 02:09:39 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1463 +1463 @@
     psF64 zSum = 1.0;
 
+    psF64 ***vsX;
+    psF64 **vsY;
+    psF64 *vsZ;
+
     nXsum = 2*nXterm;
     nYsum = 2*nYterm;
@@ -1477 +1481 @@
     // 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;
+    xSum = 1.0;
+    vsX = sums;
+    for (int i = 0; i < nXsum; i++, vsX++) {
+        ySum = xSum;
+        vsY = *vsX;
+        for (int j = 0; j < nYsum; j++, vsY++) {
+            zSum = ySum;
+            vsZ = *vsY;
+            for (int k = 0; k < nZsum; k++, vsZ++) {
+                // sums[i][j][k] = xSum;
+                *vsZ = zSum;
+                zSum *= z;
             }
             ySum *= y;
         }
-        zSum *= z;
+        xSum *= x;
     }
     return (sums);
@@ -2636 +2644 @@
     psImageInit (A, 0.0);
 
+    // use these pointers to minimize the indirections
+    psF64 *vB = B->data.F64;
+    psF64 **vA = A->data.F64;
+    psU8 ***vM = myPoly->mask;
+
     // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)*y, ...
+
+    # define NTERM(X,Y,Z)((X)*nYterm*nZterm + (Y)*nZterm + (Z))
 
     // Build the B and A data structs.
@@ -2658 +2673 @@
             for (int iy = 0; iy < nYterm; iy++) {
                 for (int iz = 0; iz < nZterm; iz++) {
-                    if (myPoly->mask[ix][iy][iz])
+                    if (vM[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;
+                    int nx = NTERM(ix,iy,iz);
+                    vB[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
                 }
             }
@@ -2669 +2684 @@
             for (int iy = 0; iy < nYterm; iy++) {
                 for (int iz = 0; iz < nZterm; iz++) {
-                    if (myPoly->mask[ix][iy][iz])
+                    if (vM[ix][iy][iz])
                         continue;
-                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    for (int jx = 0; jx < nXterm; jx++) {
-                        for (int jy = 0; jy < nYterm; jy++) {
-                            for (int jz = 0; jz < nZterm; jz++) {
-                                if (myPoly->mask[jx][jy][jz])
+                    int nx = NTERM(ix,iy,iz);
+                    for (int jx = ix; jx < nXterm; jx++) {
+                        for (int jy = (jx == ix) ? 0 : 0; jy < nYterm; jy++) {
+                            for (int jz = (jy == iy) ? 0 : 0; jz < nZterm; jz++) {
+                                if (vM[jx][jy][jz])
                                     continue;
-                                int ny = jx+jy*nXterm+jz*nXterm*nYterm;
-                                A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
+                                int ny = NTERM(jx,jy,jz);
+                                vA[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
                             }
                         }
@@ -2687 +2702 @@
     }
 
+    // fill in upper half of matrix A
+    for (int nx = 0; nx < nTerm; nx++) {
+        for (int ny = 0; ny < nx; ny++) {
+            vA[nx][ny] = vA[ny][nx];
+        }
+    }
+
     for (int ix = 0; ix < nXterm; ix++) {
         for (int iy = 0; iy < nYterm; iy++) {
             for (int iz = 0; iz < nZterm; iz++) {
-                if (!myPoly->mask[ix][iy][iz])
+                if (!vM[ix][iy][iz])
                     continue;
-                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                B->data.F64[nx] = 0;
+                int nx = NTERM(ix,iy,iz);
+                vB[nx] = 0;
                 for (int jx = 0; jx < nXterm; jx++) {
                     for (int jy = 0; jy < nYterm; jy++) {
                         for (int jz = 0; jz < nZterm; jz++) {
-                            int ny = jx+jy*nXterm+jz*nXterm*nYterm;
-                            A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
+                            int ny = NTERM(jx,jy,jz);
+                            vA[nx][ny] = (nx == ny) ? 1 : 0;
                         }
                     }
@@ -2713 +2735 @@
         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]);
+                int nx = NTERM(ix,iy,iz);
+                myPoly->coeff[ix][iy][iz] = vB[nx];
+                myPoly->coeffErr[ix][iy][iz] = sqrt(vA[nx][nx]);
             }
         }
@@ -2782 +2804 @@
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     if (x->type.type != PS_TYPE_F64) {
         PS_VECTOR_GEN_F64_FROM_F32(x64, x);
@@ -2793 +2816 @@
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     if (y->type.type != PS_TYPE_F64) {
         PS_VECTOR_GEN_F64_FROM_F32(y64, y);
@@ -2804 +2828 @@
     PS_ASSERT_VECTOR_NON_NULL(z, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(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;