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) $
8 ***********************************************************
11 DELTA = ---------------------------------------------------
12 G#*C**T*((SA*M#*A*M + SB*B**L)*C**T + SC*C**(T-1))*G
14 AOPT=1 AOPT=2 AOPT=3 AOPT=4
18 2. TWONRM=R#*R SAME SAME SAME
19 3. GMA=TWONRM/PRTN 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
31 14. DELTA= ------ SAME SAME SAME
34 TWONRM AND DELTA ONLY CG IS FOR CONSTANT
35 ARE NOT COMPUTED COMPUTED FOR DELTA, DELTA
36 FOR CONSTANT CONSTANT IS NOT
51 void quad_class::deltac(
57 REAL* delta, // corr. hstau
58 REAL* gma, // corr. hstau
65 static REAL twonrm; // set static. hstau
66 static REAL prtn;// set static. hstau
79 if(iterx == 1) twonrm = 0.0; //corr. hstau
82 if(!((bopt <= 2) && (popt <= 0))) {
89 // MULTIFLY R BY S FOR SPLITTING
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];
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];
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
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
121 // AT FIRST ITERATION FOR CONJUGATE GRADIANT SET G = R
125 for(i = 0; i < area; i++) {
131 if(bopt == 3) *gma = twonrm/prtn; //gma = twonrm/prtn;
135 for(i = 0; i < area; i++) {
136 g[i] = r[i] + (*gma) * g[i];
137 // g[i] = r[i] + gma * g[i];
145 // NO SPLITTING REQUIRED
147 for(i = 0; i < area; i++) {
154 for(i = 0; i < area; i++) {
162 adsmos(cg, t-1, cw1, cw2, cw3, cbcon);
165 badsmos(cg, t-1, cw1, cw2, cbcon);
168 // NO NEED FOR W=C**(T-1)*G IF SC=0
172 if(sc > Consts.zero) { //
176 for(i = 0; i < area; i++) {
184 adsmos(cg,1, cw1, cw2, cw3, cbcon);
187 badsmos(cg,1, cw1, cw2, cbcon);
190 // RETURN FOR CONSTANT DELTA AND AOPT=3
192 if((aopt == 3) && (bopt <= 2)) return;
196 mtamx(cg, v, list, weight);
198 if(sc > Consts.zero) {
202 for(i = 0; i < area; i++) {
209 if(sb > Consts.zero) {
211 for(i = 0; i < area; i++) {
216 adsmos(w, l, bw1, bw2, bw3, bbcon);
219 badsmos(w, l, bw1, bw2, bbcon);
224 for(i = 0; i < area; i++) {
229 for(i = 0; i < area; i++) {
237 if(!((aopt == 1) || (aopt == 3))) {
243 adsmos(v, t, cw1, cw2, cw3, cbcon);
246 badsmos(v, t, cw1, cw2, cbcon);
249 // FOR CONSTANT DELTA RETURN
251 if(bopt <= 2) return;
259 // INDIRECT NO SPLITTING
261 for(i = 0; i < area; i++) {
262 if(Blob.pix_basis || (!Blob.pix_basis && i%2 == Blob.pr)) {
270 // INDIRECT WITH SPLITTING
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];
282 // DIRECT METHOD CG = C**(T/2)*S**2*G
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];
292 // CHECK FOR ZERO DENOMINATOR
295 if(fabs(gtqg) > Consts.zero) {
299 *delta = twonrm / gtqg;
303 fprintf(output, "\n***zero denominator when computing delta iteratively, program stops***"); //not selfexplained hstau