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_prsemi.cpp $
5 $LastChangedRevision: 89 $
6 $Date: 2014-07-02 17:24:53 -0400 (Wed, 02 Jul 2014) $
8 ***********************************************************
10 THIS ROUTINE CALCULATES THE PARAMETERS FOR THE RICHARDSON METHOD
11 AND THE SEMI-ITERATIVE METHOD
14 Z = 4*((LAMBDA + ESTMIN)/(LAMBDA - ESTMIN))**2
16 ETA = 2/(LAMBDA + ESTMIN)
18 ORDER IS THE VECTOR CONTAINING THE ORDER FOR PICKING THE DELTA
20 DEL IS THE VECTOR OF DELTA FOR THE RICHARDSON METHOD
32 void quad_class::prsemi(REAL estmin, REAL lambda, BOOLEAN* alg, REAL* delta)
44 //REAL* del; //wei,3/2005 bug 174, wei, 10/2005
47 eta = lambda + estmin; //bug in Fortran. undefined eta hstau
48 if (fabs(z) <= Consts.zero)
52 "\n\n ***difference of eigenvalues is less than ZERO, program stops***");
57 if (!((period > 0) && (period <= 32)))
60 "\n\n ***period out of range (1-32), program stops***");
67 // RICHARDSON IS A STATIONARY METHOD
68 *delta = (REAL) 2.0 / eta;
70 "\n\n with period=1 this is stationary method with delta = %15.6e",
80 // LEBEDEV -FINOGENOV ORDERING
82 // CHECK TO SEE THAT PERIOD IS A POWER OF 2
85 for (i = 1; i < 6; i++)
95 "\n\n ***for lebedev-finogenov ordering period must");
96 fprintf(output, "\n be a power of 2, program stops***");
102 for (j = 1; j <= i; j++)
103 { // index j is not explicit in the loops.hstau
105 //for(m = 0; m < k; m++) {
106 for (m = 1; m <= k; m++)
109 order[2 * n + 1] = ll - order[n]; // corr.hstau
110 iorder = 2 * n; //corr. hstau
111 order[iorder] = order[n];
116 fprintf(output, "\n lebedev-finogenov ordering");
117 for (i = 0; i < period; i++)
121 fprintf(output, "\n ");
123 fprintf(output, "%3i ", order[i]);
128 // USER SUPPLIED ORDER
131 // ic = new INTEGER[period];
132 ic = new INTEGER[period + 1]; // corr. hstau
134 //for(i = 0; i < period; i++) {
135 for (i = 1; i <= period; i++)
141 for (i = 0; i < period; i++)
143 ///iorder = order + i;
146 if (!((order[i] > 0) && (order[i] <= period)))
149 "\n ***%2i th order out of range, program stops***",
155 for (i = 1; i <= period; i++)
161 "\n warning*** the %3ith delta is used %3i times",
166 delete[] ic; // bug 92 - Lajos - 03/02/2005
168 fprintf(output, "\n\n user supplied ordering");
169 for (i = 0; i < period; i++)
173 fprintf(output, "\n ");
175 fprintf(output, "%3i ", order[i]);
182 for (i = 0; i < period; i++)
185 order[i] = i + 1; //corr.hstau
188 fprintf(output, "\n\n natural ordering");
191 del = new REAL[period];
193 pdm = Consts.pi / ((REAL) 2.0 * ((REAL) (period))); //pi? hstau
195 for (i = 0; i < period; i++)
197 botm = eta - z * (REAL) cos((2.0 * ((REAL) (i)) - 1.0) * pdm);
198 if (fabs(botm) <= Consts.zero)
202 "\n\n ***cannot compute delta, zero denominator, program stops");
203 // if(del != NULL) delete[] del; //wei, 3/2005 bug 174, wei, 10/2005
206 del[i] = (REAL) 2.0 / botm;
209 //if(del != NULL) delete[] del; //wei, 3/2005 bug 174, wei, 10/2005