src/snark/quad_eigval.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_eigval.cpp $
   5  $LastChangedRevision: 80 $
   6  $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10      THIS SUBROUTINE ESTIMATES THE LOWER AND UPPER EIGEN VALUE
  11      FOR THE MATRIX
  12
  13            C**P(SA*C*M#*A*M + SB*C*B**L + SC)
  14
  15      WHERE
  16
  17          P = T - 1
  18
  19                           OR FOR THE MATRIX
  20
  21                      D*M#*A*M*D
  22
  23      Y IS A WORKING STORAGE ASSUMED TO BE THE RECONSTRUCTION AREA
  24 */
  25
  26 #include <cstdio>
  27 #include <cmath>
  28
  29 #include "blkdta.h"
  30 #include "geom.h"
  31 #include "consts.h"
  32 #include "uiod.h"
  33 #include <DIGGauss.h>
  34 #include "blob.h"
  35
  36 #include "quad.h"
  37
  38 void quad_class::eigval(
  39 REAL* y,
  40 BOOLEAN nomin,
  41 REAL* d,
  42 INTEGER* list,
  43 REAL* weight,
  44 REAL tol,
  45 REAL* estmin, //corr. hstau
  46                 REAL* lambda, // corr.hstau
  47                 REAL* s)
  48 {
  49
  50         REAL max, lower, upper;
  51
  52         INTEGER i;
  53         REAL* x;
  54         REAL presmn;
  55         INTEGER maxpt;
  56         REAL xnorml;
  57         REAL test;
  58
  59         INTEGER area;
  60         if (Blob.pix_basis)
  61         {
  62                 area = GeoPar.area;
  63         }
  64         else
  65         {
  66                 area = Blob.area;
  67         }
  68
  69         //     ALLOCATE STORAGE AREA
  70
  71         x = new REAL[area];
  72
  73 //     INITIALIZE X
  74
  75         for (i = 0; i < area; i++)
  76         {
  77                 x[i] = 1.0;
  78         }
  79
  80         /*
  81          FOR METHOD OF ESTIMATING MAX EIGENVALUE SEE
  82          A COMPUTER IMPLEMENTATION OF BAYESIAN ANALYSIS
  83          OF IMAGE RECONSTRUCTION
  84
  85          G. T. HERMAN AND A. LENT
  86
  87          INFORMATION AND CONTROL 31,364-384(1976)
  88
  89          CALCULATE Y = MATRIX*X
  90          */
  91         for (;;)
  92         {
  93                 matrix(y, x, d, list, weight, s);
  94
  95                 // SET UPPER AND LOWER BOUNDS
  96                 lower = Consts.infin;
  97                 upper = -Consts.infin;
  98                 *estmin = 0.0;
  99                 // COMPUTE LAMBDA
 100
 101                 *lambda = 0.0;
 102
 103                 for (i = 0; i < area; i++)
 104                 {
 105                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 106                         {
 107                                 (*lambda) += y[i];
 108
 109                                 // COMPUTE UPPER AND LOWER BOUND CHECK  FOR POSSIBLE DIVISION BY ZERO
 110
 111                                 if (x[i] <= Consts.zero)
 112                                 {
 113                                         fprintf(output, "\n\n      ***warning***");
 114                                         fprintf(output, "\n      z(%5i) = %15.6e", i, x[i]);
 115                                         fprintf(output,
 116                                                         "\n      lower and upper bound are ignored");
 117                                 }
 118
 119                                 else
 120                                 {
 121                                         test = y[i] / x[i];
 122
 123                                         // CHECK FOR NONPOSITIVE RESULTS
 124
 125                                         if (test <= 0.0)
 126                                         {
 127                                                 fprintf(output, "\n\n      ***warning***");
 128                                                 fprintf(output,
 129                                                                 "\n      negative value in matrix*x(%5i) = %e",
 130                                                                 i, y[i]);
 131                                                 fprintf(output, "\n      valve ignored");
 132                                         }
 133                                         else
 134                                         {
 135                                                 if (test < lower)
 136                                                         lower = test;
 137                                                 if (test > upper)
 138                                                         upper = test;
 139                                         }
 140                                 }
 141                         }
 142                 }
 143
 144                 if (Blob.pix_basis)
 145                 {
 146                         (*lambda) /= ((REAL) (area));
 147                 }
 148                 else
 149                 {
 150                         (*lambda) /= ((REAL) (area / 2.0));
 151                 }
 152
 153                 if (popt >= 1)
 154                 {
 155                         fprintf(output,
 156                                         "\n\n      lower estimate for largest eigenvalue %15.6e",
 157                                         lower);
 158                         fprintf(output,
 159                                         "\n            estimate for largest eigenvalue %15.6e",
 160                                         *lambda);
 161                         fprintf(output,
 162                                         "\n      upper estimate for largest eigenvalue %15.6e",
 163                                         upper);        //bug 147, wei, 2005/10
 164                 }
 165
 166                 // TEST FOR TERMINATION
 167
 168                 if (((upper - lower) / upper) < tol)
 169                 {
 170                         break;
 171                 }
 172
 173                 // EVALUATE NEXT X
 174
 175                 for (i = 0; i < area; i++)
 176                 {
 177                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 178                         {
 179                                 x[i] = y[i] / (*lambda);
 180                         }
 181                 }
 182         }
 183
 184         *estmin = upper;
 185
 186         if (!nomin)
 187         {
 188
 189                 /*
 190                  ESTIMATE SMALLEST EIGEN VALUE
 191                  FOR METHOD SEE ABOVE REFERENCE AND
 192
 193                  A FIRST COURSE IN NUMERICAL ANALYSIS
 194
 195                  A. RALSTON
 196
 197                  MCGRAW-HILL 1965 PAGE 475
 198
 199                  INITIALIZE X
 200                  */
 201
 202                 for (i = 0; i < area; i++)
 203                 {
 204
 205                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 206                         {
 207                                 x[i] = Gauss(0.0, 1.0);
 208                         }
 209                 }
 210
 211                 // CALCULATE Y = UPPER*X - MATRIX*X
 212
 213                 for (;;)
 214                 {
 215                         max = -Consts.infin;
 216                         matrix(y, x, d, list, weight, s);
 217                         for (i = 0; i < area; i++)
 218                         {
 219                                 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 220                                 {
 221                                         y[i] = upper * x[i] - y[i];
 222                                 }
 223                         }
 224
 225                         presmn = *estmin;
 226
 227                         if (!((t == 0) || (copt != 4)))
 228                         {
 229
 230                                 // NON SYMETRIC MATRIX
 231
 232                                 // FIND MAX OF x[i]
 233
 234                                 for (i = 0; i < area; i++)
 235                                 {
 236                                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 237                                         {
 238                                                 if (x[i] > max)
 239                                                 {
 240                                                         max = x[i];
 241                                                         maxpt = i;
 242                                                 }
 243                                         }
 244                                 }
 245
 246                                 // NEW ESTIMATES FOR MINIMUM EIGENVALUE
 247
 248                                 *estmin = y[maxpt] / max;
 249                         }
 250                         else
 251                         {
 252                                 //    SYMETRIC MATRIX
 253
 254                                 max = 0;
 255                                 *estmin = 0.0;
 256                                 for (i = 0; i < area; i++)
 257                                 {
 258                                         if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 259                                         {
 260                                                 max += x[i] * x[i];
 261                                                 (*estmin) += y[i] * x[i];
 262                                         }
 263                                 }
 264                                 (*estmin) /= max;
 265                         }
 266
 267                         if (popt >= 1)
 268                         {
 269                                 fprintf(output, "\n\n      estimate for largest eigenvalue of upper*i-matrix %15.6e",
 270                                                 *estmin); //bug 147, wei, 2005/10
 271                         }
 272
 273                         xnorml = *estmin;
 274                         if (fabs(*estmin) <= Consts.zero)
 275                                 xnorml = 1.0;
 276                         if (fabs(*estmin - presmn) / xnorml <= tol)
 277                         {
 278                                 break;
 279                         }
 280
 281                         for (i = 0; i < area; i++)
 282                         {
 283                                 if (Blob.pix_basis || (!Blob.pix_basis && i % 2 == Blob.pr))
 284                                 {
 285                                         x[i] = y[i] / xnorml;
 286                                 }
 287                         }
 288                 }
 289         }
 290         *estmin = upper - (*estmin);
 291
 292         if (popt >= 1)
 293         {
 294                 fprintf(output, "\n\n      estimate for smallest eigenvalue %15.6e", *estmin);  //bug 147, wei, 2005/10
 295         }
 296
 297         delete[] x;  // bug 92 - Lajos - 03/02/2005
 298
 299         return;
 300 }