src/snark/quad_deltac.cpp

   1 /*
   2  ***********************************************************
   3  $SNARK_Header: S N A R K  1 4 - A PICTURE RECONSTRUCTION PROGRAM $
   4  $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/quad_deltac.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10                                  R#*R
  11         DELTA = ---------------------------------------------------
  12                 G#*C**T*((SA*M#*A*M + SB*B**L)*C**T + SC*C**(T-1))*G
  13
  14     AOPT=1             AOPT=2            AOPT=3           AOPT=4
  15     (R=G=V)            (R=G)             (R=V)
  16
  17 1.                                       PRTNRM=TWONRM    SAME
  18 2.  TWONRM=R#*R        SAME              SAME             SAME
  19 3.                                       GMA=TWONRM/PRTN  SAME
  20 4.                                       G=R-GMA*G        SAME
  21 5.  W=C**(T-1)*G       SAME              SAME             SAME
  22 6.  CG=C*W             SAME              SAME             SAME
  23 7.  V=SA*M#*A*M*CG     SAME              SAME             SAME
  24 8.  V=V+SC*W           SAME              SAME             SAME
  25 9.  W=B**L*CG          SAME              SAME             SAME
  26 10. V=V+SB*W           SAME              SAME             SAME
  27 11.                    V=C**T*V                           SAME
  28 12.                    GTQG=V*G                           SAME
  29 13. GTQG=V*CG                            SAME
  30            TWONRM
  31 14. DELTA= ------      SAME              SAME             SAME
  32             GTQG
  33                        NOTE:             NOTE:            NOTE:
  34                        TWONRM AND DELTA  ONLY CG IS       FOR CONSTANT
  35                        ARE NOT COMPUTED  COMPUTED FOR     DELTA, DELTA
  36                        FOR CONSTANT      CONSTANT         IS NOT
  37                        DELTA             DELTA            COMPUTED
  38 */
  39
  40 #include <cstdio>
  41 #include <cmath>
  42
  43 #include "blkdta.h"
  44 #include "geom.h"
  45 #include "consts.h"
  46 #include "uiod.h"
  47 #include "blob.h"
  48
  49 #include "quad.h"
  50
  51 void quad_class::deltac(
  52   REAL* r,
  53   REAL* g,
  54   REAL* v,
  55   REAL* cg,
  56   REAL* w,
  57   REAL* delta,  // corr. hstau
  58   REAL* gma,    // corr. hstau
  59   INTEGER* list,
  60   REAL* weight,
  61   BOOLEAN* alg,
  62   REAL* s
  63 )
  64 {
  65   static REAL twonrm; // set static. hstau
  66   static REAL prtn;// set static. hstau
  67   REAL gtqg;
  68   INTEGER i;
  69
  70   INTEGER area;
  71
  72   if(Blob.pix_basis) {
  73     area = GeoPar.area;
  74   }
  75   else {
  76     area = Blob.area;
  77   }
  78
  79   if(iterx == 1) twonrm = 0.0; //corr. hstau
  80
  81
  82   if(!((bopt <= 2) && (popt <= 0))) {
  83     prtn = twonrm;
  84
  85 //C2.
  86
  87     twonrm = 0.0;
  88
  89     // MULTIFLY R BY S FOR SPLITTING
  90
  91     if(gopt > 0) {
  92       for(i = 0; i < area; i++) {
  93         if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
  94           twonrm += s[i] * r[i] * r[i];
  95         }
  96       }
  97     }
  98     else {
  99
 100       for(i = 0; i < area; i++) {
 101         if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
 102          twonrm += r[i] * r[i];
 103         }
 104       }
 105     }
 106
 107     if(popt >= 1) {
 108       //fprintf(output, "\n\n      norm square of gradient vector is %15.6e", twonrm);
 109          fprintf(output, "\n\n      squared norm of the residual is %15.6e", twonrm);      //wei, bug 191, 12/2005
 110     }
 111
 112     if(twonrm <= Consts.zero) {
 113       //fprintf(output, "\n\n      norm of gradient vector is zero, iteration stops");
 114      fprintf(output, "\n\n       norm of the residual is zero, iteration stops");  //wei, bug 191, 12/2005
 115
 116       *alg = TRUE;
 117       return;
 118     }
 119   }
 120
 121   // AT FIRST ITERATION FOR CONJUGATE GRADIANT SET G = R
 122
 123   if(aopt > 2) {
 124     if(iterx == 1) {
 125       for(i = 0; i < area; i++) {
 126         g[i] = r[i];
 127       }
 128     }
 129 //C3.
 130     else {
 131       if(bopt == 3) *gma = twonrm/prtn; //gma = twonrm/prtn;
 132
 133 //C4.
 134
 135       for(i = 0; i < area; i++) {
 136         g[i] = r[i] + (*gma) * g[i];
 137         // g[i] = r[i] + gma * g[i];
 138       }
 139     }
 140   }
 141   // INITILIZE CG TO G
 142
 143   if(gopt <= 0) {
 144
 145     // NO SPLITTING REQUIRED
 146
 147     for(i = 0; i < area; i++) {
 148       cg[i] = g[i];
 149     }
 150   }
 151   else {
 152     // SPLITTING
 153
 154     for(i = 0; i < area; i++) {
 155       cg[i] = g[i] * s[i];
 156     }
 157   }
 158
 159   if(t > 1) {
 160     if(copt != 5) {
 161       if(Blob.pix_basis) {
 162         adsmos(cg, t-1, cw1, cw2, cw3, cbcon);
 163       }
 164       else {
 165         badsmos(cg, t-1, cw1, cw2, cbcon);
 166       }
 167     }
 168     // NO NEED FOR W=C**(T-1)*G IF SC=0
 169   }
 170
 171
 172   if(sc > Consts.zero) { //
 173
 174 //C5.
 175
 176     for(i = 0; i < area; i++) {
 177       w[i] = cg[i];
 178     }
 179   }
 180 //C6.
 181
 182   if(copt != 5) {
 183     if(Blob.pix_basis) {
 184       adsmos(cg,1, cw1, cw2, cw3, cbcon);
 185     }
 186     else {
 187       badsmos(cg,1, cw1, cw2, cbcon);
 188     }
 189   }
 190   // RETURN FOR CONSTANT DELTA AND AOPT=3
 191
 192   if((aopt == 3) && (bopt <= 2)) return;
 193
 194 //C7.
 195
 196   mtamx(cg, v, list, weight);
 197
 198   if(sc > Consts.zero) {
 199
 200 //C8.
 201
 202     for(i = 0; i < area; i++) {
 203       v[i] += sc * w[i];
 204     }
 205   }
 206 //C9.
 207
 208
 209   if(sb > Consts.zero) {
 210     if(l > 0) {
 211       for(i = 0; i < area; i++) {
 212         w[i] = cg[i];
 213       }
 214
 215       if(Blob.pix_basis) {
 216         adsmos(w, l, bw1, bw2, bw3, bbcon);
 217       }
 218       else {
 219         badsmos(w, l, bw1, bw2, bbcon);
 220       }
 221
 222 //C10.
 223
 224       for(i = 0; i < area; i++) {
 225         v[i] += sb * w[i];
 226       }
 227     }
 228     else {
 229       for(i = 0; i < area; i++) {
 230         v[i] += sb * cg[i];
 231       }
 232     }
 233   }
 234
 235   // CHECK FOR METHOD
 236
 237   if(!((aopt == 1) || (aopt == 3))) {
 238
 239 //C11.
 240
 241     if(copt != 5) {
 242       if(Blob.pix_basis) {
 243         adsmos(v, t, cw1, cw2, cw3, cbcon);
 244       }
 245       else {
 246         badsmos(v, t, cw1, cw2, cbcon);
 247       }
 248     }
 249     // FOR CONSTANT DELTA RETURN
 250
 251     if(bopt <= 2) return;
 252
 253 //C12.
 254
 255     gtqg = 0.0;
 256
 257     if(gopt <= 0) {
 258
 259       // INDIRECT NO SPLITTING
 260
 261       for(i = 0; i < area; i++) {
 262         if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
 263         gtqg += v[i] * g[i];
 264         }
 265       }
 266
 267     }
 268     else {
 269
 270       // INDIRECT WITH SPLITTING
 271
 272       for(i = 0; i < area; i++) {
 273         if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
 274           gtqg += v[i] * g[i] * s[i];
 275         }
 276       }
 277     }
 278   }
 279   else {
 280 //C13.
 281
 282 // DIRECT METHOD CG = C**(T/2)*S**2*G
 283
 284     gtqg = 0.0;
 285
 286     for(i = 0; i < area; i++) {
 287       if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
 288         gtqg += v[i] * cg[i];
 289       }
 290     }
 291   }
 292   // CHECK FOR ZERO DENOMINATOR
 293
 294
 295   if(fabs(gtqg) > Consts.zero) {
 296
 297 //C14.
 298
 299     *delta = twonrm / gtqg;
 300
 301     return;
 302   }
 303   fprintf(output, "\n***zero denominator when computing delta iteratively, program stops***"); //not selfexplained hstau
 304   *alg = TRUE;
 305   return;
 306 }