Index: /branches/eam_rel9_b1/psLib/src/math/psMinimize.c =================================================================== --- /branches/eam_rel9_b1/psLib/src/math/psMinimize.c (revision 5808) +++ /branches/eam_rel9_b1/psLib/src/math/psMinimize.c (revision 5809) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.147.4.1.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-12-18 18:52:55 $ + * @version $Revision: 1.147.4.1.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-12-19 00:51:36 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1764,10 +1764,7 @@ B = psVectorAlloc(nTerm, PS_TYPE_F64); - // // Initialize data structures. - // XXX: Use psLib function. - // - PS_VECTOR_SET_F64(B, 0.0); - PS_IMAGE_SET_F64(A, 0.0); + psVectorInit (B, 0.0); + psImageInit (A, 0.0); // xSums look like: 1, x, x^2, ... x^(2n+1) @@ -1823,5 +1820,5 @@ // GaussJordan version - if (0) { + if (1) { // does the solution in place psGaussJordan (A, B); @@ -2130,10 +2127,7 @@ B = psVectorAlloc(nTerm, PS_TYPE_F64); - // // Initialize data structures. - // XXX: Use psLib function. - // - PS_VECTOR_SET_F64(B, 0.0); - PS_IMAGE_SET_F64(A, 0.0); + psVectorInit (B, 0.0); + psImageInit (A, 0.0); // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1) @@ -2168,9 +2162,11 @@ if (myPoly->mask[i][j]) continue; + int nx = i+j*nXterm; 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; + int ny = n+m*nXterm; + A->data.F64[nx][ny] += Sums->data.F64[i+n][j+m] * wt; } } @@ -2183,10 +2179,10 @@ if (!myPoly->mask[i][j]) continue; - ix = i+j*nXterm; - B->data.F64[ix] = 0; + int nx = i+j*nXterm; + B->data.F64[nx] = 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; + int ny = n+m*nXterm; + A->data.F64[nx][ny] = (nx == ny) ? 1 : 0; } } @@ -2664,5 +2660,6 @@ 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; + int nx = ix+iy*nXterm+iz*nXterm*nYterm; + B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt; } } @@ -2674,4 +2671,5 @@ if (myPoly->mask[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++) { @@ -2679,5 +2677,6 @@ 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; + int ny = jx+jy*nXterm+jz*nXterm*nYterm; + A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt; } } @@ -2693,10 +2692,10 @@ if (!myPoly->mask[ix][iy][iz]) continue; - nx = ix+iy*nXterm+iz*nXterm*nYterm; + int 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; + int ny = jx+jy*nXterm+jz*nXterm*nYterm; A->data.F64[nx][ny] = (nx == ny) ? 1 : 0; } @@ -2714,6 +2713,8 @@ for (int iy = 0; iy < nYterm; iy++) { 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]); + int nx = ix+iy*nXterm+iz*nXterm*nYterm; + int ny = jx+jy*nXterm+jz*nXterm*nYterm; + myPoly->coeff[ix][iy][iz] = B->data.F64[nx]; + myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][ny]); } } @@ -3111,5 +3112,6 @@ 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; + int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; + B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt; } } @@ -3123,4 +3125,5 @@ if (myPoly->mask[ix][iy][iz][it]) continue; + int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; for (int jx = 0; jx < nXterm; jx++) { for (int jy = 0; jy < nYterm; jy++) { @@ -3129,6 +3132,6 @@ 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; + int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm; + A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt; } } @@ -3147,5 +3150,5 @@ if (!myPoly->mask[ix][iy][iz][it]) continue; - nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; + int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; B->data.F64[nx] = 0; for (int jx = 0; jx < nXterm; jx++) { @@ -3153,5 +3156,5 @@ 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; + int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm; A->data.F64[nx][ny] = (nx == ny) ? 1 : 0; } @@ -3172,5 +3175,8 @@ for (int iz = 0; iz < nZterm; iz++) { for (int it = 0; it < nTterm; it++) { - myPoly->coeff[ix][iy][iz][it] = B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+iz*nXterm*nYterm*nZterm]; + int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm; + int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm; + myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx]; + myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][ny]); } }