src/snark/unkps.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/unkps.cpp $
   5  $LastChangedRevision: 122 $
   6  $Date: 2014-07-09 18:18:59 -0400 (Wed, 09 Jul 2014) $
   7  $Author: agulati $
   8  ***********************************************************
   9
  10  */
  11
  12 #include <cstdio>
  13 #include <cmath>
  14
  15 #include "blkdta.h"
  16
  17 #include "unkps.h"
  18
  19
  20 void unkps( INTEGER imax, INTEGER istep, REAL* xs, REAL* x, INTEGER m, INTEGER nipf, INTEGER krad, INTEGER mc, INTEGER ir, INTEGER n, INTEGER* ix, INTEGER* ipu, INTEGER nfixp, INTEGER k4m4, BOOLEAN virt)
  21 {
  22         BOOLEAN dovirt;
  23
  24         REAL md4;
  25         REAL kradm;
  26         INTEGER k4mx;
  27         INTEGER k4mxm;
  28         INTEGER k4minc;
  29         INTEGER jlim;
  30         INTEGER i;
  31         INTEGER j;
  32         INTEGER js;
  33         INTEGER inc;
  34         INTEGER jl;
  35         INTEGER nn;
  36         INTEGER kp;
  37         INTEGER kp1;
  38         INTEGER ju;
  39         INTEGER jus;
  40         INTEGER kx;
  41         INTEGER ia;
  42         INTEGER l;
  43         INTEGER ind;
  44         INTEGER inx;
  45         INTEGER ib;
  46         INTEGER ind2;
  47         INTEGER kumx;
  48         INTEGER ku;
  49         INTEGER kvmn;
  50         INTEGER kv;
  51         INTEGER kp2;
  52         INTEGER k4mxm1;
  53         INTEGER kp1x;
  54         INTEGER kp2x;
  55         INTEGER kp1xx;
  56         INTEGER kp2xx;
  57         INTEGER kpl;
  58         INTEGER kount;
  59
  60         dovirt = virt && (istep > 1);
  61         md4 = m / krad;
  62         if ((md4 == 0) && (nipf < 2))
  63                 return;
  64         kradm = krad - 1;
  65         k4mx = k4m4 * ir;
  66         k4mxm = k4mx - k4m4;
  67         k4minc = k4m4;
  68         ix[0] = 0;
  69         jlim = 1;
  70
  71         for (i = 0; i < md4; i++)
  72         {
  73                 for (j = 1; j <= jlim; j++)
  74                 {
  75                         js = j - 1;
  76                         inc = 0;
  77                         ix[j - 1] *= krad;
  78                         for (jl = 0; jl < kradm; jl++)
  79                         {
  80                                 js += jlim;
  81                                 inc += k4mx;
  82                                 ix[js] = ix[j - 1] + inc;
  83                         }
  84                 }
  85                 jlim *= krad;
  86         }
  87
  88         nn = n / (1 << m);
  89
  90         if (!(nipf < 2))
  91         {
  92                 k4minc = k4mxm;
  93
  94                 for (kp = 0; kp < k4m4; kp++)
  95                 {
  96                         kp1 = ix[kp] + kp;
  97                         ju = nfixp;
  98
  99                         do
 100                         {
 101                                 jus = ju;
 102                                 kx = kp1 + ipu[ju];
 103                                 ia = kx * istep;
 104
 105
 106                                 for (l = 0; l < imax; l++)
 107                                 {
 108                                         ind = l + ia;
 109                                         CA(xs, 0, l)= CA(x, 0, ind);
 110                                         CA(xs, 1, l)= CA(x, 1, ind);
 111                                 }
 112
 113                                 ju++;
 114
 115                                 do
 116                                 {
 117                                         inx = kp1 + ipu[ju];
 118                                         ia = kx * istep;
 119
 120                                         ib = inx * istep;
 121
 122                                         for (l = 0; l < imax; l++)
 123                                         {
 124                                                 ind = l + ia;
 125                                                 ind2 = l + ib;
 126                                                 CA(x, 0, ind)= CA(x, 0, ind2);
 127                                                 CA(x, 1, ind)= CA(x, 1, ind2);
 128                                         }
 129
 130                                         ju++;
 131
 132                                         kx = inx;
 133
 134                                 } while (ipu[ju] != 0);
 135
 136                                 ia = kx * istep;
 137
 138                                 for (l = 0; l < imax; l++)
 139                                 {
 140                                         ind = l + ia;
 141                                         CA(x, 0, ind)= CA(xs, 0, l);
 142                                         CA(x, 1, ind)= CA(xs, 1, l);
 143                                 }
 144
 145                                 ju++;
 146
 147                         } while (ipu[ju] < n);
 148                 }
 149         }
 150
 151         if (md4 < 1)
 152                 return;
 153
 154         kumx = k4m4 - 1;
 155
 156         for (ku = 0; ku < kumx; ku++)
 157         {
 158                 kvmn = ku + 1;
 159                 for (kv = kvmn; kv < k4m4; kv++)
 160                 {
 161
 162                         kp1 = ix[ku] + kv;
 163                         kp2 = ix[kv] + ku;
 164
 165                         k4mxm1 = k4mxm + 1;
 166
 167                         for (inc = 0; inc < k4mxm1; inc += k4minc)
 168                         {
 169
 170                                 kp1x = kp1 + inc;
 171                                 kp2x = kp2 + inc;
 172
 173                                 ia = kp1x * istep;
 174
 175                                 ib = kp2x * istep;
 176
 177
 178                                 for (l = 0; l < imax; l++)
 179                                 {
 180                                         kp1xx = l + ia;
 181                                         kp2xx = l + ib;
 182
 183                                         CA(xs, 0, l)= CA(x, 0, kp1xx);
 184                                         CA(xs, 1, l)= CA(x, 1, kp1xx);
 185
 186                                         CA(x, 0, kp1xx)= CA(x, 0, kp2xx);
 187                                         CA(x, 1, kp1xx)= CA(x, 1, kp2xx);
 188
 189                                         CA(x, 0, kp2xx)= CA(xs, 0, l);
 190                                         CA(x, 1, kp2xx)= CA(xs, 1, l);
 191                                 }
 192                         }
 193
 194                         // code below is not tested !!! mk
 195                         // after some correction it seems to work - Laszlo
 196
 197                         if (!(nipf < 2))
 198                         {
 199
 200                                 if (!(nfixp == 0))
 201                                 {
 202
 203                                         for (ju = 0; ju < nfixp; ju++)
 204                                         {
 205
 206                                                 ind = kp1 + ipu[ju];
 207                                                 ia = (ind) * istep;
 208
 209
 210                                                 ind = kp2 + ipu[ju];
 211                                                 ib = (ind) * istep;
 212
 213
 214                                                 for (l = 0; l < imax; l++)
 215                                                 {
 216                                                         kp1x = l + ia;
 217                                                         kp2x = l + ib;
 218
 219                                                         CA(xs, 0, l)= CA(x, 0, kp1x);
 220                                                         CA(xs, 1, l)= CA(x, 1, kp1x);
 221
 222
 223                                                         CA(x, 0, kp1x)= CA(x, 0, kp2x);
 224                                                         CA(x, 1, kp1x)= CA(x, 1, kp2x);
 225
 226                                                         CA(x, 0, kp2x)= CA(xs, 0, l);
 227                                                         CA(x, 1, kp2x)= CA(xs, 1, l);
 228                                                 }
 229                                         }
 230                                 }
 231
 232                                 ju = nfixp;
 233                                 do
 234                                 {
 235                                         kount = 1;
 236                                         jus = ju;
 237                                         kpl = 0;
 238                                         kx = kp1 + ipu[ju];
 239
 240                                         for(;;)
 241                                         {
 242                                                 ia = (kx) * istep;
 243
 244                                                 for(l = 0; l < imax; l++)
 245                                                 {
 246                                                         ind = l + ia;
 247                                                         CA(xs, 0, l) = CA(x, 0, ind);
 248                                                         CA(xs, 1, l) = CA(x, 1, ind);
 249                                                 }
 250
 251                                                 ju++;
 252
 253                                                 do
 254                                                 {
 255                                                         do
 256                                                         {
 257                                                                 kpl = 1 - kpl;
 258
 259                                                                 if(kpl == 1) inx = kp2 + ipu[ju];
 260                                                                 if(kpl == 0) inx = kp1 + ipu[ju];
 261
 262                                                                 ia = (kx) * istep;
 263
 264
 265                                                                 ib = (inx) * istep;
 266
 267
 268                                                                 for(l = 0; l < imax; l++)
 269                                                                 {
 270
 271                                                                         ind = l + ia;
 272                                                                         ind2 = l + ib;
 273
 274                                                                         CA(x, 0, ind) = CA(x, 0, ind2);
 275                                                                         CA(x, 1, ind) = CA(x, 1, ind2);
 276
 277                                                                 }
 278
 279                                                                 ju++;
 280                                                                 kx = inx;
 281
 282                                                         }while(ipu[ju] != 0);
 283
 284                                                         if(kount == 2) goto L12;
 285
 286                                                         kount = 2;
 287                                                         ju = jus;
 288
 289                                                 }while(kpl == 0);
 290
 291                                                 ia = (kx) * istep;
 292
 293
 294                                                 for(l = 0; l < imax; l++)
 295                                                 {
 296                                                         ind = l + ia;
 297
 298                                                         CA(x, 0, ind) = CA(xs, 0, l);
 299                                                         CA(x, 1, ind) = CA(xs, 1, l);
 300                                                 }
 301
 302                                                 kx = kp2 + ipu[ju];
 303                                         }
 304
 305                                         L12:
 306
 307                                         ia = (kx) * istep;
 308
 309                                         for(l = 0; l < imax; l++)
 310                                         {
 311                                                 ind = l + ia;
 312
 313                                                 CA(x, 0, ind) = CA(xs, 0, l);
 314                                                 CA(x, 1, ind) = CA(xs, 1, l);
 315                                         }
 316
 317                                         ju++;
 318
 319                                 }while(ipu[ju] < n);
 320                         }
 321                 }
 322         }
 323         return;
 324 }