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_eigval.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 THIS SUBROUTINE ESTIMATES THE LOWER AND UPPER EIGEN VALUE
13 C**P(SA*C*M#*A*M + SB*C*B**L + SC)
23 Y IS A WORKING STORAGE ASSUMED TO BE THE RECONSTRUCTION AREA
38 void quad_class::eigval(
45 REAL* estmin, //corr. hstau
46 REAL* lambda, // corr.hstau
50 REAL max, lower, upper;
69 // ALLOCATE STORAGE AREA
75 for (i = 0; i < area; i++)
81 FOR METHOD OF ESTIMATING MAX EIGENVALUE SEE
82 A COMPUTER IMPLEMENTATION OF BAYESIAN ANALYSIS
83 OF IMAGE RECONSTRUCTION
85 G. T. HERMAN AND A. LENT
87 INFORMATION AND CONTROL 31,364-384(1976)
89 CALCULATE Y = MATRIX*X
93 matrix(y, x, d, list, weight, s);
95 // SET UPPER AND LOWER BOUNDS
97 upper = -Consts.infin;
103 for (i = 0; i < area; i++)
105 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
109 // COMPUTE UPPER AND LOWER BOUND CHECK FOR POSSIBLE DIVISION BY ZERO
111 if (x[i] <= Consts.zero)
113 fprintf(output, "\n\n ***warning***");
114 fprintf(output, "\n z(%5i) = %15.6e", i, x[i]);
116 "\n lower and upper bound are ignored");
123 // CHECK FOR NONPOSITIVE RESULTS
127 fprintf(output, "\n\n ***warning***");
129 "\n negative value in matrix*x(%5i) = %e",
131 fprintf(output, "\n valve ignored");
146 (*lambda) /= ((REAL) (area));
150 (*lambda) /= ((REAL) (area / 2.0));
156 "\n\n lower estimate for largest eigenvalue %15.6e",
159 "\n estimate for largest eigenvalue %15.6e",
162 "\n upper estimate for largest eigenvalue %15.6e",
163 upper); //bug 147, wei, 2005/10
166 // TEST FOR TERMINATION
168 if (((upper - lower) / upper) < tol)
175 for (i = 0; i < area; i++)
177 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
179 x[i] = y[i] / (*lambda);
190 ESTIMATE SMALLEST EIGEN VALUE
191 FOR METHOD SEE ABOVE REFERENCE AND
193 A FIRST COURSE IN NUMERICAL ANALYSIS
197 MCGRAW-HILL 1965 PAGE 475
202 for (i = 0; i < area; i++)
205 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
207 x[i] = Gauss(0.0, 1.0);
211 // CALCULATE Y = UPPER*X - MATRIX*X
216 matrix(y, x, d, list, weight, s);
217 for (i = 0; i < area; i++)
219 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
221 y[i] = upper * x[i] - y[i];
227 if (!((t == 0) || (copt != 4)))
230 // NON SYMETRIC MATRIX
234 for (i = 0; i < area; i++)
236 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
246 // NEW ESTIMATES FOR MINIMUM EIGENVALUE
248 *estmin = y[maxpt] / max;
256 for (i = 0; i < area; i++)
258 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
261 (*estmin) += y[i] * x[i];
269 fprintf(output, "\n\n estimate for largest eigenvalue of upper*i-matrix %15.6e",
270 *estmin); //bug 147, wei, 2005/10
274 if (fabs(*estmin) <= Consts.zero)
276 if (fabs(*estmin - presmn) / xnorml <= tol)
281 for (i = 0; i < area; i++)
283 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
285 x[i] = y[i] / xnorml;
290 *estmin = upper - (*estmin);
294 fprintf(output, "\n\n estimate for smallest eigenvalue %15.6e", *estmin); //bug 147, wei, 2005/10
297 delete[] x; // bug 92 - Lajos - 03/02/2005