IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 18, 2005, 2:51:36 PM (21 years ago)
Author:
magnier
Message:

finished up coeffErr and coeffMask code

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5808 r5809  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.147.4.1.2.2 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-12-18 18:52:55 $
     12 *  @version $Revision: 1.147.4.1.2.3 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-12-19 00:51:36 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    17641764    B     = psVectorAlloc(nTerm, PS_TYPE_F64);
    17651765
    1766     //
    17671766    // Initialize data structures.
    1768     // XXX: Use psLib function.
    1769     //
    1770     PS_VECTOR_SET_F64(B, 0.0);
    1771     PS_IMAGE_SET_F64(A, 0.0);
     1767    psVectorInit (B, 0.0);
     1768    psImageInit (A, 0.0);
    17721769
    17731770    // xSums look like: 1, x, x^2, ... x^(2n+1)
    … …  
    18231820
    18241821    // GaussJordan version
    1825     if (0) {
     1822    if (1) {
    18261823        // does the solution in place
    18271824        psGaussJordan (A, B);
    … …  
    21302127    B = psVectorAlloc(nTerm, PS_TYPE_F64);
    21312128
    2132     //
    21332129    // Initialize data structures.
    2134     // XXX: Use psLib function.
    2135     //
    2136     PS_VECTOR_SET_F64(B, 0.0);
    2137     PS_IMAGE_SET_F64(A, 0.0);
     2130    psVectorInit (B, 0.0);
     2131    psImageInit (A, 0.0);
    21382132
    21392133    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
    … …  
    21682162                if (myPoly->mask[i][j])
    21692163                    continue;
     2164                int nx = i+j*nXterm;
    21702165                for (int n = 0; n < nXterm; n++) {
    21712166                    for (int m = 0; m < nYterm; m++) {
    21722167                        if (myPoly->mask[n][m])
    21732168                            continue;
    2174                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
     2169                        int ny = n+m*nXterm;
     2170                        A->data.F64[nx][ny] += Sums->data.F64[i+n][j+m] * wt;
    21752171                    }
    21762172                }
    … …  
    21832179            if (!myPoly->mask[i][j])
    21842180                continue;
    2185             ix = i+j*nXterm;
    2186             B->data.F64[ix] = 0;
     2181            int nx = i+j*nXterm;
     2182            B->data.F64[nx] = 0;
    21872183            for (int n = 0; n < nXterm; n++) {
    21882184                for (int m = 0; m < nYterm; m++) {
    2189                     iy = n+m*nXterm;
    2190                     A->data.F64[ix][iy] = (ix == iy) ? 1 : 0;
     2185                    int ny = n+m*nXterm;
     2186                    A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
    21912187                }
    21922188            }
    … …  
    26642660                    if (myPoly->mask[ix][iy][iz])
    26652661                        continue;
    2666                     B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
     2662                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
     2663                    B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz] * wt;
    26672664                }
    26682665            }
    … …  
    26742671                    if (myPoly->mask[ix][iy][iz])
    26752672                        continue;
     2673                    int nx = ix+iy*nXterm+iz*nXterm*nYterm;
    26762674                    for (int jx = 0; jx < nXterm; jx++) {
    26772675                        for (int jy = 0; jy < nYterm; jy++) {
    … …  
    26792677                                if (myPoly->mask[jx][jy][jz])
    26802678                                    continue;
    2681                                 A->data.F64[ix+iy*nXterm+iz*nXterm*nYterm][ix+iy*nXterm+iz*nXterm*nYterm] += Sums[ix+jx][iy+jy][iz+jz] * wt;
     2679                                int ny = jx+jy*nXterm+jz*nXterm*nYterm;
     2680                                A->data.F64[nx][ny] += Sums[ix+jx][iy+jy][iz+jz] * wt;
    26822681                            }
    26832682                        }
    … …  
    26932692                if (!myPoly->mask[ix][iy][iz])
    26942693                    continue;
    2695                 nx = ix+iy*nXterm+iz*nXterm*nYterm;
     2694                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
    26962695                B->data.F64[nx] = 0;
    26972696                for (int jx = 0; jx < nXterm; jx++) {
    26982697                    for (int jy = 0; jy < nYterm; jy++) {
    26992698                        for (int jz = 0; jz < nZterm; jz++) {
    2700                             ny = jx+jy*nXterm+jz*nXterm*nYterm;
     2699                            int ny = jx+jy*nXterm+jz*nXterm*nYterm;
    27012700                            A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
    27022701                        }
    … …  
    27142713        for (int iy = 0; iy < nYterm; iy++) {
    27152714            for (int iz = 0; iz < nZterm; iz++) {
    2716                 myPoly->coeff[ix][iy][iz] = B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm];
    2717                 myPoly->coeffErr[ix][iy][iz] = sqrt(B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm][ix+iy*nXterm+iz*nXterm*nYterm]);
     2715                int nx = ix+iy*nXterm+iz*nXterm*nYterm;
     2716                int ny = jx+jy*nXterm+jz*nXterm*nYterm;
     2717                myPoly->coeff[ix][iy][iz] = B->data.F64[nx];
     2718                myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][ny]);
    27182719            }
    27192720        }
    … …  
    31113112                        if (myPoly->mask[ix][iy][iz][it])
    31123113                            continue;
    3113                         B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
     3114                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
     3115                        B->data.F64[nx] += f->data.F64[k] * Sums[ix][iy][iz][it] * wt;
    31143116                    }
    31153117                }
    … …  
    31233125                        if (myPoly->mask[ix][iy][iz][it])
    31243126                            continue;
     3127                        int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    31253128                        for (int jx = 0; jx < nXterm; jx++) {
    31263129                            for (int jy = 0; jy < nYterm; jy++) {
    … …  
    31293132                                        if (myPoly->mask[jx][jy][jz][jt])
    31303133                                            continue;
    3131                                         A->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm][jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm]
    3132                                         += Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
     3134                                        int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
     3135                                        A->data.F64[nx][ny]+= Sums[ix+jx][iy+jy][iz+jz][it+jt] * wt;
    31333136                                    }
    31343137                                }
    … …  
    31473150                    if (!myPoly->mask[ix][iy][iz][it])
    31483151                        continue;
    3149                     nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
     3152                    int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
    31503153                    B->data.F64[nx] = 0;
    31513154                    for (int jx = 0; jx < nXterm; jx++) {
    … …  
    31533156                            for (int jz = 0; jz < nZterm; jz++) {
    31543157                                for (int jt = 0; jt < nTterm; jt++) {
    3155                                     ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
     3158                                    int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
    31563159                                    A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
    31573160                                }
    … …  
    31723175            for (int iz = 0; iz < nZterm; iz++) {
    31733176                for (int it = 0; it < nTterm; it++) {
    3174                     myPoly->coeff[ix][iy][iz][it] = B->data.F64[ix+iy*nXterm+iz*nXterm*nYterm+iz*nXterm*nYterm*nZterm];
     3177                    int nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
     3178                    int ny = jx+jy*nXterm+jz*nXterm*nYterm+jt*nXterm*nYterm*nZterm;
     3179                    myPoly->coeff[ix][iy][iz][it] = B->data.F64[nx];
     3180                    myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][ny]);
    31753181                }
    31763182            }
Note: See TracChangeset for help on using the changeset viewer.