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_deltad.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
11 DELTA = ---------------------
14 NOTE: R = M#*A*(M*X-P)
16 AOPT=1 AOPT=2 AOPT=3 AOPT=4
21 2. TWONRM=D**2*R**2 SAME SAME SAME
22 3. GMA=TWONRM/PRTN SAME
24 5. CG=D**2*G SAME SAME SAME
25 6. V=M#*A*M*CG SAME SAME SAME
26 7. GTQG=CG*V SAME SAME SAME
28 8. DELTA=------ SAME SAME SAME
32 TWONRM AND FOR CONSTANT
33 DELTA ARE NOT DELTA ONLY SAME
34 COMPUTED FOR V IS COMPUTED
49 void quad_class::deltad(
55 REAL* delta, // corr. hstau
56 REAL* gma, //corr. hstau
61 static REAL twonrm; // set static. hstau
62 static REAL prtn; // set static. hstau
78 twonrm = 0.0; //corr? hstau
82 if (!((bopt <= 2) && (popt <= 0)))
89 for (i = 0; i < area; i++)
91 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
93 twonrm += d[i] * r[i] * r[i];
99 fprintf(output, "\n\n squared norm of the residual is %15.6e",
100 twonrm); //wei, bug 191, 12/2005
103 if (twonrm <= Consts.zero)
106 "\n\n norm of the residual is zero, iteration stops"); //wei, bug 191, 12/2005
118 for (i = 0; i < area; i++)
126 *gma = twonrm / prtn; //gma = twonrm/prtn;
130 for (i = 0; i < area; i++)
132 g[i] = r[i] + (*gma) * g[i]; // corr. hstau
139 for (i = 0; i < area; i++)
146 // RETURN FOR AOPT .LE. 2 AND BOPT .LE. 2
148 if ((bopt <= 2) && (aopt <= 2))
149 return; // corr.hstau
152 mtamx(cg, v, list, weight);
160 for (i = 0; i < area; i++)
162 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
164 gtqg += cg[i] * v[i];
169 // CHECK FOR ZERO DENOMINATOR
171 if (fabs(gtqg) > Consts.zero)
175 *delta = twonrm / gtqg;
180 fprintf(output, "\n***zero denominator when computing delta iteratively, program stops***"); //not selfexplained hstau