IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 1, 2005, 9:55:16 AM (21 years ago)
Author:
eugene
Message:

adding ensemble photometry

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psSparse.c

    r5653 r5654  
    11# include "psphot.h"
    2 
    3 typedef struct {
    4     psVector *Aij;
    5     psVector *Bfj;
    6     psVector *Qii;
    7     psVector *Si;
    8     psVector *Sj;
    9     int Nelem;
    10     int Nrows;
    11 } psSparse;
    122
    133void psSparseMatrixTest () {
    … …  
    166    psSparse *sparse = psSparseAlloc (3, 9);
    177
    18     psSparseMatrixElement (sparse, 0, 0, 1.0);
    19     psSparseMatrixElement (sparse, 1, 1, 1.0);
     8    psSparseMatrixElement (sparse, 0, 0, 3.0);
     9    psSparseMatrixElement (sparse, 1, 1, 2.0);
    2010    psSparseMatrixElement (sparse, 2, 2, 1.0);
    2111
    22     psSparseMatrixElement (sparse, 0, 1, 0.1);
    23     psSparseMatrixElement (sparse, 0, 2, 0.1);
    24 
    25     psSparseMatrixResort (sparse);
     12    psSparseMatrixElement (sparse, 1, 0, 0.1);
     13    psSparseMatrixElement (sparse, 2, 0, -0.1);
     14
     15    psSparseResort (sparse);
    2616    for (int i = 0; i < sparse->Nelem; i++) {
    2717        fprintf (stderr, "%d %d %f\n",
    … …  
    3323    psVector *x = psVectorAlloc (3, PS_DATA_F32);
    3424    x->data.F32[0] = 3;
    35     x->data.F32[0] = 5;
    36     x->data.F32[0] = 7;
    37 
    38     psVector B = psSparseMatrixTimesVector (NULL, sparse, x);
     25    x->data.F32[1] = 5;
     26    x->data.F32[2] = 7;
     27    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     28
     29    psVector *B = psSparseMatrixTimesVector (NULL, sparse, x);
    3930    fprintf (stderr, "B: %f %f %f\n", B->data.F32[0], B->data.F32[1], B->data.F32[2]);
    4031
    … …  
    4334    sparse->Bfj->data.F32[2] = B->data.F32[2];
    4435
    45     x = psSparseSolve (x, sparse);
     36    x = psSparseSolve (x, sparse, 0);
     37    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     38
     39    x = psSparseSolve (x, sparse, 1);
     40    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     41
     42    x = psSparseSolve (x, sparse, 2);
     43    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     44
     45    x = psSparseSolve (x, sparse, 3);
     46    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
     47
     48    x = psSparseSolve (x, sparse, 4);
    4649    fprintf (stderr, "x: %f %f %f\n", x->data.F32[0], x->data.F32[1], x->data.F32[2]);
    4750    return;
    … …  
    5053void psSparseResort (psSparse *sparse) {
    5154
    52     Nelem = sparse->Nelem;
     55    int Nelem = sparse->Nelem;
    5356
    5457    psVector *index = psVectorSortIndex (NULL, sparse->Sj);
    55     psVector Aij = sparse->Aij;
    56     psVector Si = sparse->Si;
    57     psVector Sj = sparse->Sj;
     58    psVector *Aij = sparse->Aij;
     59    psVector *Si = sparse->Si;
     60    psVector *Sj = sparse->Sj;
    5861
    5962    // allocate new temporary vectors
    … …  
    6164    psVector *tSi  = psVectorAlloc (Nelem, PS_DATA_S32);
    6265    psVector *tSj  = psVectorAlloc (Nelem, PS_DATA_S32);
    63     for (i = 0; i < Nelem; i++) {
    64         j = index->data.U32[i];
     66    for (int i = 0; i < Nelem; i++) {
     67        int j = index->data.U32[i];
    6568        tAij->data.F32[i] = Aij->data.F32[j];
    6669        tSi->data.S32[i]  = Si->data.S32[j];
    … …  
    8083void psSparseMatrixElement (psSparse *sparse, int i, int j, float value) {
    8184
     85    int k;
     86
    8287    if (i < j) {
    8388        fprintf (stderr, "*** error: subdiagonal element ***\n");
    … …  
    9196        // check vectors lengths and extend if needed
    9297        if (sparse->Nelem >= sparse->Aij->nalloc) {
    93             sparse->Aij->nalloc += 100;     
    94             psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
    95             sparse->Si->nalloc += 100;     
    96             psVectorRealloc (sparse->Si, sparse->Si->nalloc);
    97             sparse->Sj->nalloc += 100;     
    98             psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
     98            psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
     99            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
     100            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
    99101        }
    100102
    … …  
    111113        // check vectors lengths and extend if needed
    112114        if (sparse->Nelem >= sparse->Aij->nalloc - 1) {
    113             sparse->Aij->nalloc += 100;     
    114             psVectorRealloc (sparse->Aij, sparse->Aij->nalloc);
    115             sparse->Si->nalloc += 100;     
    116             psVectorRealloc (sparse->Si, sparse->Si->nalloc);
    117             sparse->Sj->nalloc += 100;     
    118             psVectorRealloc (sparse->Sj, sparse->Sj->nalloc);
     115            psVectorRealloc (sparse->Aij, sparse->Aij->nalloc + 100);
     116            psVectorRealloc (sparse->Si,  sparse->Si->nalloc + 100);
     117            psVectorRealloc (sparse->Sj,  sparse->Sj->nalloc + 100);
    119118        }
    120119
    … …  
    139138void psSparseVectorElement (psSparse *sparse, int i, float value) {
    140139
    141     sparse->Bjf->data.F32[i] = value;
     140    sparse->Bfj->data.F32[i] = value;
    142141    return;
    143142}
    … …  
    146145psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector) {
    147146
     147    int i, Nelem;
     148    float F;
     149
    148150    if (output == NULL) {
    149         psVector *output = psVectorAlloc (vector->n, PS_DATA_F32);
     151        output = psVectorAlloc (vector->n, PS_DATA_F32);
    150152    }
    151153
    … …  
    154156        F = 0;
    155157        while (matrix->Sj->data.S32[Nelem] == j) {
    156             i = matrix->Sj->data.S32[Nelem];
    157             F += vector->data.F32[i] * matrix->Aij->data.F32[i];
     158            i = matrix->Si->data.S32[Nelem];
     159            F += vector->data.F32[i] * matrix->Aij->data.F32[Nelem];
     160            Nelem++;
    158161        }
    159162        output->data.F32[j] = F;
    … …  
    162165}
    163166
    164 psVector *psSparseSolve (psVector *guess, psSparse *sparse) {
     167psVector *psSparseSolve (psVector *guess, psSparse *sparse, int Niter) {
    165168
    166169    psF32 dG;
    … …  
    169172    psVector *Bfj = sparse->Bfj;
    170173
    171     guess = psVectorCopy (guess, Bfj);
     174    guess = psVectorCopy (guess, Bfj, PS_DATA_F32);
    172175
    173176    // temporary storage for intermediate results
    174177    psVector *dQ = psVectorAlloc (guess->n, PS_DATA_F32);
    175178
    176     for (int j = 0; j < 2; j++) {
     179    for (int j = 0; j < Niter; j++) {
    177180        dQ = psSparseMatrixTimesVector (dQ, sparse, guess);
    178         for (int i = 0; i < dG->n; i++) {
     181        for (int i = 0; i < dQ->n; i++) {
    179182            dG = (dQ->data.F32[i] - Bfj->data.F32[i]) / Qii->data.F32[i];
    180183            guess->data.F32[i] -= dG;
    … …  
    187190static void psSparseFree (psSparse *sparse) {
    188191    if (sparse == NULL) return;
    189     psFree (sparse.Aij);
    190     psFree (sparse.Bfj);
    191     psFree (sparse.Qii);
    192     psFree (sparse.Si);
    193     psFree (sparse.Sj);
     192    psFree (sparse->Aij);
     193    psFree (sparse->Bfj);
     194    psFree (sparse->Qii);
     195    psFree (sparse->Si);
     196    psFree (sparse->Sj);
    194197    return;
    195198}
    … …  
    212215
    213216    sparse->Nrows = Nrows;
     217    psMemSetDeallocator(sparse, (psFreeFunc) psSparseFree);
    214218    return (sparse);
    215219}
Note: See TracChangeset for help on using the changeset viewer.