src/snark/quad_prsemi.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_prsemi.cpp $
   5  $LastChangedRevision: 89 $
   6  $Date: 2014-07-02 17:24:53 -0400 (Wed, 02 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  THIS ROUTINE CALCULATES THE PARAMETERS FOR THE RICHARDSON METHOD
  11  AND THE SEMI-ITERATIVE METHOD
  12  HERE
  13
  14  Z = 4*((LAMBDA + ESTMIN)/(LAMBDA - ESTMIN))**2
  15
  16  ETA = 2/(LAMBDA + ESTMIN)
  17
  18  ORDER IS THE VECTOR CONTAINING THE ORDER FOR PICKING THE DELTA
  19  FOR RICHARDSON METHOD
  20  DEL IS THE VECTOR OF DELTA FOR THE RICHARDSON METHOD
  21  */
  22
  23 #include <cstdio>
  24 #include <cmath>
  25
  26 #include "blkdta.h"
  27 #include "consts.h"
  28 #include "uiod.h"
  29
  30 #include "quad.h"
  31
  32 void quad_class::prsemi(REAL estmin, REAL lambda, BOOLEAN* alg, REAL* delta)
  33 {
  34         INTEGER i;
  35         INTEGER j;
  36         INTEGER k;
  37         INTEGER n;
  38         INTEGER m;
  39         INTEGER ll;
  40         INTEGER* ic;
  41         INTEGER iorder;
  42         REAL pdm;
  43         REAL botm;
  44         //REAL* del;   //wei,3/2005  bug 174, wei, 10/2005
  45
  46         z = lambda - estmin;
  47         eta = lambda + estmin;  //bug in Fortran. undefined eta hstau
  48         if (fabs(z) <= Consts.zero)
  49         {
  50                 *alg = TRUE;
  51                 fprintf(output,
  52                                 "\n\n       ***difference of eigenvalues is less than ZERO, program stops***");
  53                 return;
  54         }
  55
  56         // RICHARDSON METHOD
  57         if (!((period > 0) && (period <= 32)))
  58         {
  59                 fprintf(output,
  60                                 "\n\n       ***period out of range (1-32), program stops***");
  61                 *alg = TRUE;
  62                 return;
  63         }
  64
  65         if (period == 1)
  66         {
  67                 // RICHARDSON IS A STATIONARY METHOD
  68                 *delta = (REAL) 2.0 / eta;
  69                 fprintf(output,
  70                                 "\n\n       with period=1 this is stationary method with delta = %15.6e",
  71                                 *delta);
  72                 bopt = 1;
  73                 return;
  74         }
  75
  76         if ((dopt - 2) >= 0)
  77         {
  78                 if ((dopt - 2) > 0)
  79                 {
  80                         // LEBEDEV -FINOGENOV ORDERING
  81
  82                         // CHECK TO SEE THAT PERIOD IS A POWER OF 2
  83
  84                         m = period;
  85                         for (i = 1; i < 6; i++)
  86                         {   //corr. hstau
  87                                 if (m != (m / 2) * 2)
  88                                         goto L5;
  89                                 m = m / 2;
  90                                 if (m == 1)
  91                                         goto L6;
  92                         }
  93
  94                         L5: fprintf(output,
  95                                         "\n\n      ***for lebedev-finogenov ordering period must");
  96                         fprintf(output, "\n      be a power of 2, program stops***");
  97                         *alg = TRUE;
  98                         return;
  99
 100                         L6: order[0] = 1;
 101                         k = 1;
 102                         for (j = 1; j <= i; j++)
 103                         {  // index j is not explicit in the loops.hstau
 104                                 ll = 2 * k + 1;
 105                                 //for(m = 0; m < k; m++) {
 106                                 for (m = 1; m <= k; m++)
 107                                 {
 108                                         n = k - m;
 109                                         order[2 * n + 1] = ll - order[n];   // corr.hstau
 110                                         iorder = 2 * n;  //corr. hstau
 111                                         order[iorder] = order[n];
 112                                 }
 113                                 k = 2 * k;
 114                         }
 115
 116                         fprintf(output, "\n      lebedev-finogenov ordering");
 117                         for (i = 0; i < period; i++)
 118                         {
 119                                 if ((i % 32) == 0)
 120                                 {
 121                                         fprintf(output, "\n      ");
 122                                 }
 123                                 fprintf(output, "%3i ", order[i]);
 124                         }
 125                 }
 126                 else
 127                 {
 128                         // USER SUPPLIED ORDER
 129                         // CHECK USER ORDER
 130
 131                         // ic = new INTEGER[period];
 132                         ic = new INTEGER[period + 1]; // corr. hstau
 133
 134                         //for(i = 0; i < period; i++) {
 135                         for (i = 1; i <= period; i++)
 136                         {  // corr.hstau
 137                                 ///ici = ic + i;
 138                                 ic[i] = 0;
 139                         }
 140
 141                         for (i = 0; i < period; i++)
 142                         {
 143                                 ///iorder = order + i;
 144                                 j = order[i];
 145                                 ic[j]++;
 146                                 if (!((order[i] > 0) && (order[i] <= period)))
 147                                 {
 148                                         fprintf(output,
 149                                                         "\n      ***%2i th order out of range, program stops***",
 150                                                         i);
 151                                         *alg = TRUE;
 152                                 }
 153                         }
 154
 155                         for (i = 1; i <= period; i++)
 156                         {  // corr.hstau
 157                                 ///ici = ic + i;
 158                                 if (ic[i] != 1)
 159                                 {
 160                                         fprintf(output,
 161                                                         "\n       warning*** the %3ith delta is used %3i times",
 162                                                         i, ic[i]);
 163                                 }
 164                         }
 165
 166                         delete[] ic;  // bug 92 - Lajos - 03/02/2005
 167
 168                         fprintf(output, "\n\n      user supplied ordering");
 169                         for (i = 0; i < period; i++)
 170                         {
 171                                 if ((i % 32) == 0)
 172                                 {
 173                                         fprintf(output, "\n      ");
 174                                 }
 175                                 fprintf(output, "%3i ", order[i]);
 176                         }
 177                 }
 178         }
 179         // NATURAL ORDER
 180         else
 181         {
 182                 for (i = 0; i < period; i++)
 183                 {
 184                         //order[i] = i;
 185                         order[i] = i + 1; //corr.hstau
 186                 }
 187
 188                 fprintf(output, "\n\n      natural ordering");
 189                 // COMPUTE DELTA
 190         }
 191         del = new REAL[period];
 192
 193         pdm = Consts.pi / ((REAL) 2.0 * ((REAL) (period))); //pi? hstau
 194
 195         for (i = 0; i < period; i++)
 196         {
 197                 botm = eta - z * (REAL) cos((2.0 * ((REAL) (i)) - 1.0) * pdm);
 198                 if (fabs(botm) <= Consts.zero)
 199                 {
 200                         *alg = TRUE;
 201                         fprintf(output,
 202                                         "\n\n      ***cannot compute delta, zero denominator, program stops");
 203                         // if(del != NULL) delete[] del;  //wei, 3/2005  bug 174, wei, 10/2005
 204                         return;
 205                 }
 206                 del[i] = (REAL) 2.0 / botm;
 207
 208         }
 209         //if(del != NULL) delete[] del;  //wei, 3/2005 bug 174, wei, 10/2005
 210         return;
 211 }