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

Timestamp:

Dec 18, 2005, 8:52:55 AM (21 years ago)

Author:

magnier

Message:

added masked orders and coeff errors

File:

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

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

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.147.4.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 20:26:29 $
+ *  @version $Revision: 1.147.4.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-18 18:52:55 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1775 +1775 @@
     // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects
     // XXX EAM : this function is only valid for data vectors of F64
+
     for (int k = 0; k < f->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] && maskValue)) {
@@ -1792 +1793 @@
         }
         for (int i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
 
         // we could skip half of the array and assign at the end
-        // we must handle masked orders
+        // skip the masked orders, they get assigned (0 or 1) below
         for (int i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             for (int j = 0; j < nTerm; j++) {
+                if (myPoly->mask[j])
+                    continue;
                 A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
             }
         }
     }
+
+    // set masked elements in A,B appropriately
+    for (int i = 0; i < nTerm; i++) {
+        if (!myPoly->mask[i])
+            continue;
+        B->data.F64[i] = 0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = (i == j) ? 1 : 0;
+        }
+    }
+
 
     // GaussJordan version
@@ -1813 +1831 @@
         for (int k = 0; k < nTerm; k++) {
             myPoly->coeff[k] = B->data.F64[k];
+            myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
         }
     } else {
@@ -1825 +1844 @@
         for (int k = 0; k < nTerm; k++) {
             myPoly->coeff[k] = coeffs->data.F64[k];
+            // XXX EAM : LUD does not provide inverse(A)
+            // myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
         }
         psFree(ALUD);
@@ -2137 +2158 @@
         for (int n = 0; n < nXterm; n++) {
             for (int m = 0; m < nYterm; m++) {
+                if (myPoly->mask[n][m])
+                    continue;
                 B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
@@ -2143 +2166 @@
         for (int i = 0; i < nXterm; i++) {
             for (int j = 0; j < nYterm; j++) {
+                if (myPoly->mask[i][j])
+                    continue;
                 for (int n = 0; n < nXterm; n++) {
                     for (int m = 0; m < nYterm; m++) {
+                        if (myPoly->mask[n][m])
+                            continue;
                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
                     }
+                }
+            }
+        }
+    }
+
+    for (int i = 0; i < nXterm; i++) {
+        for (int j = 0; j < nYterm; j++) {
+            if (!myPoly->mask[i][j])
+                continue;
+            ix = i+j*nXterm;
+            B->data.F64[ix] = 0;
+            for (int n = 0; n < nXterm; n++) {
+                for (int m = 0; m < nYterm; m++) {
+                    iy = n+m*nXterm;
+                    A->data.F64[ix][iy] = (ix == iy) ? 1 : 0;
                 }
             }
@@ -2159 +2201 @@
         for (int m = 0; m < nYterm; m++) {
             myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+            myPoly->coeffErr[n][m] = sqrt(A->data.F64[n+m*nXterm][n+m*nXterm]);
         }
     }
@@ -2619 +2662 @@
             for (int iy = 0; iy < nYterm; iy++) {
                 for (int iz = 0; iz < nZterm; iz++) {
+                    if (myPoly->mask[ix][iy][iz])
+                        continue;
                     B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
                 }
@@ -2627 +2672 @@
             for (int iy = 0; iy < nYterm; iy++) {
                 for (int iz = 0; iz < nZterm; iz++) {
+                    if (myPoly->mask[ix][iy][iz])
+                        continue;
                     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])
+                                    continue;
                                 A->data.F64[ix+iy*nXterm+iz*nXterm*nYterm][ix+iy*nXterm+iz*nXterm*nYterm] += Sums[ix+jx][iy+jy][iz+jz] * wt;
                             }
@@ -2639 +2688 @@
     }
 
+    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])
+                    continue;
+                nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                B->data.F64[nx] = 0;
+                for (int jx = 0; jx < nXterm; jx++) {
+                    for (int jy = 0; jy < nYterm; jy++) {
+                        for (int jz = 0; jz < nZterm; jz++) {
+                            ny = jx+jy*nXterm+jz*nXterm*nYterm;
+                            A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
     // does the solution in place
     psGaussJordan (A, B);
@@ -2647 +2715 @@
             for (int iz = 0; iz < nZterm; iz++) {
                 myPoly->coeff[ix][iy][iz] = B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm];
+                myPoly->coeffErr[ix][iy][iz] = sqrt(B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm][ix+iy*nXterm+iz*nXterm*nYterm]);
             }
         }
@@ -3040 +3109 @@
                 for (int iz = 0; iz < nZterm; iz++) {
                     for (int it = 0; it < nTterm; it++) {
+                        if (myPoly->mask[ix][iy][iz][it])
+                            continue;
                         B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
                     }
@@ -3050 +3121 @@
                 for (int iz = 0; iz < nZterm; iz++) {
                     for (int it = 0; it < nTterm; it++) {
+                        if (myPoly->mask[ix][iy][iz][it])
+                            continue;
                         for (int jx = 0; jx < nXterm; jx++) {
                             for (int jy = 0; jy < nYterm; jy++) {
                                 for (int jz = 0; jz < nZterm; jz++) {
                                     for (int jt = 0; jt < nTterm; jt++) {
-                                        A->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm][ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm]
+                                        if (myPoly->mask[jx][jy][jz][jt])
+                                            continue;
+                                        A->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm][jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm]
                                         += Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
                                     }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    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++) {
+                    if (!myPoly->mask[ix][iy][iz][it])
+                        continue;
+                    nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                    B->data.F64[nx] = 0;
+                    for (int jx = 0; jx < nXterm; jx++) {
+                        for (int jy = 0; jy < nYterm; jy++) {
+                            for (int jz = 0; jz < nZterm; jz++) {
+                                for (int jt = 0; jt < nTterm; jt++) {
+                                    ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
+                                    A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                                 }
                             }