Index: /branches/eam_rel9_b1/psLib/src/math/psMinimize.c =================================================================== --- /branches/eam_rel9_b1/psLib/src/math/psMinimize.c (revision 5807) +++ /branches/eam_rel9_b1/psLib/src/math/psMinimize.c (revision 5808) @@ -10,6 +10,6 @@ * @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,4 +1775,5 @@ // 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,15 +1793,32 @@ } 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,4 +1831,5 @@ 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,4 +1844,6 @@ 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,4 +2158,6 @@ 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,8 +2166,27 @@ 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,4 +2201,5 @@ 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,4 +2662,6 @@ 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,7 +2672,11 @@ 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,4 +2688,23 @@ } + 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,4 +2715,5 @@ 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,4 +3109,6 @@ 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,11 +3121,38 @@ 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; } }