Robust/src/Benchmarks/Prefetch/2DFFT/dsm/fft2d.java

   1 public class fft2d extends Thread {
   2   //Title:        2-d mixed radix FFT.
   3   //Version:
   4   //Copyright:    Copyright (c) 1998
   5   //Author:       Dongyan Wang
   6   //Company:      University of Wisconsin-Milwaukee.
   7   //Description:
   8   //              . Use fft1d to perform fft2d.
   9   //
  10   // Code borrowed from :Java Digital Signal Processing book by Lyon and Rao
  11
  12   public Matrix data1, data2;
  13   public int x0, x1;
  14
  15   // Constructor: 2-d FFT of Complex data.
  16   public fft2d(Matrix data1, Matrix data2, int x0, int x1) {
  17     this.data1 = data1;
  18     this.data2 = data2;
  19     this.x0 = x0;
  20     this.x1 = x1;
  21   }
  22
  23   public void run() {
  24     fft1d fft1, fft2;
  25     Barrier barr;
  26     barr = new Barrier("128.195.175.84");
  27     double tempdataRe[][];
  28     double tempdataIm[][];
  29     int rowlength, columnlength;
  30     int start, end;
  31
  32     // Calculate FFT for each row of the data.
  33     atomic {
  34       rowlength = data1.M;
  35       columnlength = data1.N;
  36       tempdataRe = data1.dataRe;
  37       tempdataIm = data1.dataIm;
  38       start = x0;
  39       end = x1;
  40       fft1 = new fft1d(columnlength);
  41       fft2 = new fft1d(rowlength);
  42       for (int i = x0; i < x1; i++) {
  43         //input of FFT
  44         double inputRe[] = tempdataRe[i]; //local array
  45         double inputIm[] = tempdataIm[i];
  46         fft(fft1, inputRe, inputIm);
  47       } //end of for
  48     }
  49
  50     //Start Barrier
  51     Barrier.enterBarrier(barr);
  52
  53     // Tranpose data.
  54     if (start == 0) {
  55       atomic {
  56         transpose(tempdataRe,tempdataIm, data2.dataRe,data2.dataIm, rowlength, columnlength);
  57       }
  58     }
  59
  60     //Start Barrier
  61     Barrier.enterBarrier(barr);
  62
  63     // Calculate FFT for each column of the data.
  64     double transtempRe[][];
  65     double transtempIm[][];
  66     atomic {
  67       transtempRe = data2.dataRe;
  68       transtempIm = data2.dataIm;
  69       for (int j = start; j < end; j++) {
  70         //input of FFT
  71         double inputRe[] = transtempRe[j]; //local array
  72         double inputIm[] = transtempIm[j];
  73         fft(fft2, inputRe, inputIm);
  74       } //end of fft2 for
  75     }
  76   } //end of run
  77
  78   public void transpose(double[][] tempdataRe, double[][] tempdataIm, double[][] outputRe,
  79                         double[][] outputIm, int rowlength, int columnlength) {
  80     for(int i = 0; i<rowlength; i++) {
  81       double tRe[] = tempdataRe[i];
  82       double tIm[] = tempdataIm[i];
  83       for(int j = 0; j<columnlength; j++) {
  84         outputRe[j][i] = tRe[j];
  85         outputIm[j][i] = tIm[j];
  86       }
  87     }
  88   }
  89
  90   public static void main(String[] args) {
  91     int NUM_THREADS = 1;
  92     int SIZE = 800;
  93     int inputWidth = 10;
  94     if(args.length>0) {
  95       NUM_THREADS=Integer.parseInt(args[0]);
  96       if(args.length > 1)
  97         SIZE = Integer.parseInt(args[1]);
  98     }
  99
 100     System.printString("Num threads = " + NUM_THREADS + " SIZE= " + SIZE + "\n");
 101
 102     int[] mid = new int[8];
 103     mid[0] = (128<<24)|(195<<16)|(175<<8)|84; //dw-10
 104     mid[1] = (128<<24)|(195<<16)|(175<<8)|85; //dw-11
 105     mid[2] = (128<<24)|(195<<16)|(175<<8)|86; //dw-12
 106     mid[3] = (128<<24)|(195<<16)|(175<<8)|87; //dw-13
 107     mid[4] = (128<<24)|(195<<16)|(175<<8)|88; //dw-14
 108     mid[5] = (128<<24)|(195<<16)|(175<<8)|89; //dw-15
 109     mid[6] = (128<<24)|(195<<16)|(175<<8)|90; //dw-16
 110     mid[7] = (128<<24)|(195<<16)|(175<<8)|91; //dw-17
 111
 112     // Start Barrier Server
 113     BarrierServer mybarr;
 114     atomic {
 115       mybarr = global new BarrierServer(NUM_THREADS);
 116     }
 117     mybarr.start(mid[0]);
 118
 119     Matrix data1;
 120     Matrix data2;
 121
 122     // Create threads to do FFT
 123     fft2d[] myfft2d;
 124     atomic {
 125       // Set up data for FFT transform
 126       data1 = global new Matrix(SIZE, SIZE);
 127       data2 = global new Matrix(SIZE, SIZE);
 128       data1.setValues(); //Input Matrix
 129       data2.setZeros(); //Transpose Matrix
 130       myfft2d = global new fft2d[NUM_THREADS];
 131       int increment = SIZE/NUM_THREADS;
 132       int base = 0;
 133       for(int i =0 ; i<NUM_THREADS; i++) {
 134         if((i+1)==NUM_THREADS)
 135           myfft2d[i] = global new fft2d(data1, data2, base, SIZE);
 136         else
 137           myfft2d[i] = global new fft2d(data1, data2, base, base+increment);
 138         base+=increment;
 139       }
 140     }
 141
 142     boolean waitfordone=true;
 143     while(waitfordone) {
 144       atomic {
 145         if (mybarr.done)
 146           waitfordone=false;
 147       }
 148     }
 149
 150     fft2d tmp;
 151     //Start a thread to compute each c[l,n]
 152     for(int i = 0; i<NUM_THREADS; i++) {
 153       atomic {
 154         tmp = myfft2d[i];
 155       }
 156       tmp.start(mid[i]);
 157     }
 158
 159     //Wait for thread to finish
 160     for(int i = 0; i<NUM_THREADS; i++) {
 161       atomic {
 162         tmp = myfft2d[i];
 163       }
 164       tmp.join();
 165     }
 166
 167     System.printString("2DFFT done! \n");
 168   }
 169
 170   public static void fft(fft1d myfft, double inputRe[], double inputIm[]) {
 171     //output of FFT
 172     double outputRe[] = myfft.outputRe;
 173     double outputIm[] = myfft.outputIm;
 174     // intermediate results
 175     double temRe[] = myfft.temRe;
 176     double temIm[] = myfft.temIm;
 177     //Permute() operation
 178     permute(myfft, outputRe, outputIm, inputRe, inputIm);
 179
 180     //System.printString("ready to twiddle");
 181     for (int factorIndex = 0; factorIndex < myfft.NumofFactors; factorIndex++)
 182       twiddle(factorIndex, myfft, temRe, temIm, outputRe, outputIm);
 183
 184     //System.printString("ready to copy");
 185     // Copy the output[] data to input[], so the output can be
 186     // returned in the input array.
 187     for (int i = 0; i < myfft.N; i++) {
 188       inputRe[i] = outputRe[i];
 189       inputIm[i] = outputIm[i];
 190     }
 191   }
 192
 193   private static void permute(fft1d myfft, double[] outputRe, double[] outputIm, double[] inputRe, double[] inputIm) {
 194     int count[] = new int[myfft.MaxFactorsNumber];
 195     int j;
 196     int k = 0;
 197
 198     for (int i = 0; i < myfft.N - 1; i++) {
 199       outputRe[i] = inputRe[k];
 200       outputIm[i] = inputIm[k];
 201       j = 0;
 202       k = k + myfft.remain[j];
 203       count[0] = count[0] + 1;
 204       while (count[j] >= myfft.factors[j]) {
 205         count[j] = 0;
 206         int tmp;
 207         if(j == 0)
 208           tmp = myfft.N;
 209         else
 210           tmp = myfft.remain[j - 1];
 211         k = k - tmp + myfft.remain[j + 1];
 212         j++;
 213         count[j] = count[j] + 1;
 214       }
 215     }
 216     outputRe[myfft.N - 1] = inputRe[myfft.N - 1];
 217     outputIm[myfft.N - 1] = inputIm[myfft.N - 1];
 218   }   // End of function permute().
 219
 220   private static void twiddle(int factorIndex, fft1d myfft, double[] temRe, double[] temIm,
 221                               double[] outputRe, double[] outputIm) {
 222     // Get factor data.
 223     int sofarRadix = myfft.sofar[factorIndex];
 224     int radix = myfft.factors[factorIndex];
 225     int remainRadix = myfft.remain[factorIndex];
 226
 227     double tem;   // Temporary variable to do data exchange.
 228
 229     double W = 2 * (double) Math.setPI() / (sofarRadix * radix);
 230     double cosW = (double) Math.cos(W);
 231     double sinW = -(double) Math.sin(W);
 232
 233     double twiddleRe[] = new double[radix];
 234     double twiddleIm[] = new double[radix];
 235     double twRe = 1.0f, twIm = 0f;
 236
 237     //Initialize twiddle addBk.address variables.
 238     int dataOffset = 0, groupOffset = 0, address = 0;
 239
 240     for (int dataNo = 0; dataNo < sofarRadix; dataNo++) {
 241       //System.printString("datano="+dataNo);
 242       if (sofarRadix > 1) {
 243         twiddleRe[0] = 1.0f;
 244         twiddleIm[0] = 0.0f;
 245         twiddleRe[1] = twRe;
 246         twiddleIm[1] = twIm;
 247         for (int i = 2; i < radix; i++) {
 248           twiddleRe[i] = twRe * twiddleRe[i - 1] - twIm * twiddleIm[i - 1];
 249           twiddleIm[i] = twIm * twiddleRe[i - 1] + twRe * twiddleIm[i - 1];
 250         }
 251         tem = cosW * twRe - sinW * twIm;
 252         twIm = sinW * twRe + cosW * twIm;
 253         twRe = tem;
 254       }
 255       for (int groupNo = 0; groupNo < remainRadix; groupNo++) {
 256         //System.printString("groupNo="+groupNo);
 257         if ((sofarRadix > 1) && (dataNo > 0)) {
 258           temRe[0] = outputRe[address];
 259           temIm[0] = outputIm[address];
 260           int blockIndex = 1;
 261           do {
 262             address = address + sofarRadix;
 263             temRe[blockIndex] = twiddleRe[blockIndex] * outputRe[address] -
 264                                 twiddleIm[blockIndex] * outputIm[address];
 265             temIm[blockIndex] = twiddleRe[blockIndex] * outputIm[address] +
 266                                 twiddleIm[blockIndex] * outputRe[address];
 267             blockIndex++;
 268           } while (blockIndex < radix);
 269         } else {
 270           for (int i = 0; i < radix; i++) {
 271             //System.printString("temRe.length="+temRe.length);
 272             //System.printString("i = "+i);
 273             temRe[i] = outputRe[address];
 274             temIm[i] = outputIm[address];
 275             address += sofarRadix;
 276           }
 277         }
 278         //System.printString("radix="+radix);
 279         if(radix == 2) {
 280           tem = temRe[0] + temRe[1];
 281           temRe[1] = temRe[0] - temRe[1];
 282           temRe[0] = tem;
 283           tem = temIm[0] + temIm[1];
 284           temIm[1] = temIm[0] - temIm[1];
 285           temIm[0] = tem;
 286         } else if( radix == 3) {
 287           double t1Re = temRe[1] + temRe[2];
 288           double t1Im = temIm[1] + temIm[2];
 289           temRe[0] = temRe[0] + t1Re;
 290           temIm[0] = temIm[0] + t1Im;
 291
 292           double m1Re = myfft.cos2to3PI * t1Re;
 293           double m1Im = myfft.cos2to3PI * t1Im;
 294           double m2Re = myfft.sin2to3PI * (temIm[1] - temIm[2]);
 295           double m2Im = myfft.sin2to3PI * (temRe[2] - temRe[1]);
 296           double s1Re = temRe[0] + m1Re;
 297           double s1Im = temIm[0] + m1Im;
 298
 299           temRe[1] = s1Re + m2Re;
 300           temIm[1] = s1Im + m2Im;
 301           temRe[2] = s1Re - m2Re;
 302           temIm[2] = s1Im - m2Im;
 303         } else if(radix == 4) {
 304           fft4(temRe, temIm);
 305         } else if(radix == 5) {
 306           fft5(myfft, temRe, temIm);
 307         } else if(radix == 8) {
 308           fft8(myfft, temRe, temIm);
 309         } else if(radix == 10) {
 310           fft10(myfft, temRe, temIm);
 311         } else {
 312           fftPrime(radix, temRe, temIm);
 313         }
 314         address = groupOffset;
 315         for (int i = 0; i < radix; i++) {
 316           outputRe[address] = temRe[i];
 317           outputIm[address] = temIm[i];
 318           address += sofarRadix;
 319         }
 320         groupOffset += sofarRadix * radix;
 321         address = groupOffset;
 322       }
 323       groupOffset = ++dataOffset;
 324       address = groupOffset;
 325     }
 326   } //twiddle operation
 327
 328   // The two arguments dataRe[], dataIm[] are mainly for using in fft8();
 329   private static void fft4(double dataRe[], double dataIm[]) {
 330     double t1Re,t1Im, t2Re,t2Im;
 331     double m2Re,m2Im, m3Re,m3Im;
 332
 333     t1Re = dataRe[0] + dataRe[2];
 334     t1Im = dataIm[0] + dataIm[2];
 335     t2Re = dataRe[1] + dataRe[3];
 336     t2Im = dataIm[1] + dataIm[3];
 337
 338     m2Re = dataRe[0] - dataRe[2];
 339     m2Im = dataIm[0] - dataIm[2];
 340     m3Re = dataIm[1] - dataIm[3];
 341     m3Im = dataRe[3] - dataRe[1];
 342
 343     dataRe[0] = t1Re + t2Re;
 344     dataIm[0] = t1Im + t2Im;
 345     dataRe[2] = t1Re - t2Re;
 346     dataIm[2] = t1Im - t2Im;
 347     dataRe[1] = m2Re + m3Re;
 348     dataIm[1] = m2Im + m3Im;
 349     dataRe[3] = m2Re - m3Re;
 350     dataIm[3] = m2Im - m3Im;
 351   }   // End of function fft4().
 352
 353   // The two arguments dataRe[], dataIm[] are mainly for using in fft10();
 354   private static void fft5(fft1d myfft, double dataRe[], double dataIm[]) {
 355     double t1Re,t1Im, t2Re,t2Im, t3Re,t3Im, t4Re,t4Im, t5Re,t5Im;
 356     double m1Re,m1Im, m2Re,m2Im, m3Re,m3Im, m4Re,m4Im, m5Re,m5Im;
 357     double s1Re,s1Im, s2Re,s2Im, s3Re,s3Im, s4Re,s4Im, s5Re,s5Im;
 358
 359     t1Re = dataRe[1] + dataRe[4];
 360     t1Im = dataIm[1] + dataIm[4];
 361     t2Re = dataRe[2] + dataRe[3];
 362     t2Im = dataIm[2] + dataIm[3];
 363     t3Re = dataRe[1] - dataRe[4];
 364     t3Im = dataIm[1] - dataIm[4];
 365     t4Re = dataRe[3] - dataRe[2];
 366     t4Im = dataIm[3] - dataIm[2];
 367     t5Re = t1Re + t2Re;
 368     t5Im = t1Im + t2Im;
 369
 370     dataRe[0] = dataRe[0] + t5Re;
 371     dataIm[0] = dataIm[0] + t5Im;
 372
 373     m1Re = myfft.c51 * t5Re;
 374     m1Im = myfft.c51 * t5Im;
 375     m2Re = myfft.c52 * (t1Re - t2Re);
 376     m2Im = myfft.c52 * (t1Im - t2Im);
 377     m3Re = -(myfft.c53) * (t3Im + t4Im);
 378     m3Im = myfft.c53 * (t3Re + t4Re);
 379     m4Re = -(myfft.c54) * t4Im;
 380     m4Im = myfft.c54 * t4Re;
 381     m5Re = -(myfft.c55) * t3Im;
 382     m5Im = myfft.c55 * t3Re;
 383
 384     s3Re = m3Re - m4Re;
 385     s3Im = m3Im - m4Im;
 386     s5Re = m3Re + m5Re;
 387     s5Im = m3Im + m5Im;
 388     s1Re = dataRe[0] + m1Re;
 389     s1Im = dataIm[0] + m1Im;
 390     s2Re = s1Re + m2Re;
 391     s2Im = s1Im + m2Im;
 392     s4Re = s1Re - m2Re;
 393     s4Im = s1Im - m2Im;
 394
 395     dataRe[1] = s2Re + s3Re;
 396     dataIm[1] = s2Im + s3Im;
 397     dataRe[2] = s4Re + s5Re;
 398     dataIm[2] = s4Im + s5Im;
 399     dataRe[3] = s4Re - s5Re;
 400     dataIm[3] = s4Im - s5Im;
 401     dataRe[4] = s2Re - s3Re;
 402     dataIm[4] = s2Im - s3Im;
 403   }   // End of function fft5().
 404
 405   private static void fft8(fft1d myfft, double[] temRe, double[] temIm) {
 406     double data1Re[] = new double[4];
 407     double data1Im[] = new double[4];
 408     double data2Re[] = new double[4];
 409     double data2Im[] = new double[4];
 410     double tem;
 411
 412     // To improve the speed, use direct assaignment instead for loop here.
 413     data1Re[0] = temRe[0];
 414     data2Re[0] = temRe[1];
 415     data1Re[1] = temRe[2];
 416     data2Re[1] = temRe[3];
 417     data1Re[2] = temRe[4];
 418     data2Re[2] = temRe[5];
 419     data1Re[3] = temRe[6];
 420     data2Re[3] = temRe[7];
 421
 422     data1Im[0] = temIm[0];
 423     data2Im[0] = temIm[1];
 424     data1Im[1] = temIm[2];
 425     data2Im[1] = temIm[3];
 426     data1Im[2] = temIm[4];
 427     data2Im[2] = temIm[5];
 428     data1Im[3] = temIm[6];
 429     data2Im[3] = temIm[7];
 430
 431     fft4(data1Re, data1Im);
 432     fft4(data2Re, data2Im);
 433
 434     tem = myfft.OnetoSqrt2 * (data2Re[1] + data2Im[1]);
 435     data2Im[1] = myfft.OnetoSqrt2 * (data2Im[1] - data2Re[1]);
 436     data2Re[1] = tem;
 437     tem = data2Im[2];
 438     data2Im[2] = -data2Re[2];
 439     data2Re[2] = tem;
 440     tem = myfft.OnetoSqrt2 * (data2Im[3] - data2Re[3]);
 441     data2Im[3] = -(myfft.OnetoSqrt2) * (data2Re[3] + data2Im[3]);
 442     data2Re[3] = tem;
 443
 444     temRe[0] = data1Re[0] + data2Re[0];
 445     temRe[4] = data1Re[0] - data2Re[0];
 446     temRe[1] = data1Re[1] + data2Re[1];
 447     temRe[5] = data1Re[1] - data2Re[1];
 448     temRe[2] = data1Re[2] + data2Re[2];
 449     temRe[6] = data1Re[2] - data2Re[2];
 450     temRe[3] = data1Re[3] + data2Re[3];
 451     temRe[7] = data1Re[3] - data2Re[3];
 452
 453     temIm[0] = data1Im[0] + data2Im[0];
 454     temIm[4] = data1Im[0] - data2Im[0];
 455     temIm[1] = data1Im[1] + data2Im[1];
 456     temIm[5] = data1Im[1] - data2Im[1];
 457     temIm[2] = data1Im[2] + data2Im[2];
 458     temIm[6] = data1Im[2] - data2Im[2];
 459     temIm[3] = data1Im[3] + data2Im[3];
 460     temIm[7] = data1Im[3] - data2Im[3];
 461   }   // End of function fft8().
 462
 463   private static void fft10(fft1d myfft, double[] temRe, double[] temIm) {
 464     double data1Re[] = new double[5];
 465     double data1Im[] = new double[5];
 466     double data2Re[] = new double[5];
 467     double data2Im[] = new double[5];
 468
 469     // To improve the speed, use direct assaignment instead for loop here.
 470     data1Re[0] = temRe[0];
 471     data2Re[0] = temRe[5];
 472     data1Re[1] = temRe[2];
 473     data2Re[1] = temRe[7];
 474     data1Re[2] = temRe[4];
 475     data2Re[2] = temRe[9];
 476     data1Re[3] = temRe[6];
 477     data2Re[3] = temRe[1];
 478     data1Re[4] = temRe[8];
 479     data2Re[4] = temRe[3];
 480
 481     data1Im[0] = temIm[0];
 482     data2Im[0] = temIm[5];
 483     data1Im[1] = temIm[2];
 484     data2Im[1] = temIm[7];
 485     data1Im[2] = temIm[4];
 486     data2Im[2] = temIm[9];
 487     data1Im[3] = temIm[6];
 488     data2Im[3] = temIm[1];
 489     data1Im[4] = temIm[8];
 490     data2Im[4] = temIm[3];
 491
 492     fft5(myfft, data1Re, data1Im);
 493     fft5(myfft, data2Re, data2Im);
 494
 495     temRe[0] = data1Re[0] + data2Re[0];
 496     temRe[5] = data1Re[0] - data2Re[0];
 497     temRe[6] = data1Re[1] + data2Re[1];
 498     temRe[1] = data1Re[1] - data2Re[1];
 499     temRe[2] = data1Re[2] + data2Re[2];
 500     temRe[7] = data1Re[2] - data2Re[2];
 501     temRe[8] = data1Re[3] + data2Re[3];
 502     temRe[3] = data1Re[3] - data2Re[3];
 503     temRe[4] = data1Re[4] + data2Re[4];
 504     temRe[9] = data1Re[4] - data2Re[4];
 505
 506     temIm[0] = data1Im[0] + data2Im[0];
 507     temIm[5] = data1Im[0] - data2Im[0];
 508     temIm[6] = data1Im[1] + data2Im[1];
 509     temIm[1] = data1Im[1] - data2Im[1];
 510     temIm[2] = data1Im[2] + data2Im[2];
 511     temIm[7] = data1Im[2] - data2Im[2];
 512     temIm[8] = data1Im[3] + data2Im[3];
 513     temIm[3] = data1Im[3] - data2Im[3];
 514     temIm[4] = data1Im[4] + data2Im[4];
 515     temIm[9] = data1Im[4] - data2Im[4];
 516   }   // End of function fft10().
 517
 518   private static void fftPrime(int radix, double[] temRe, double[] temIm) {
 519     // Initial WRe, WIm.
 520     double W = 2 * (double) Math.setPI() / radix;
 521     double cosW = (double) Math.cos(W);
 522     double sinW = -(double) Math.sin(W);
 523     double WRe[] = new double[radix];
 524     double WIm[] = new double[radix];
 525
 526     WRe[0] = 1;
 527     WIm[0] = 0;
 528     WRe[1] = cosW;
 529     WIm[1] = sinW;
 530
 531     for (int i = 2; i < radix; i++) {
 532       WRe[i] = cosW * WRe[i - 1] - sinW * WIm[i - 1];
 533       WIm[i] = sinW * WRe[i - 1] + cosW * WIm[i - 1];
 534     }
 535
 536     // FFT of prime length data, using DFT, can be improved in the future.
 537     double rere, reim, imre, imim;
 538     int j, k;
 539     int max = (radix + 1) / 2;
 540
 541     double tem1Re[] = new double[max];
 542     double tem1Im[] = new double[max];
 543     double tem2Re[] = new double[max];
 544     double tem2Im[] = new double[max];
 545
 546     for (j = 1; j < max; j++) {
 547       tem1Re[j] = temRe[j] + temRe[radix - j];
 548       tem1Im[j] = temIm[j] - temIm[radix - j];
 549       tem2Re[j] = temRe[j] - temRe[radix - j];
 550       tem2Im[j] = temIm[j] + temIm[radix - j];
 551     }
 552
 553     for (j = 1; j < max; j++) {
 554       temRe[j] = temRe[0];
 555       temIm[j] = temIm[0];
 556       temRe[radix - j] = temRe[0];
 557       temIm[radix - j] = temIm[0];
 558       k = j;
 559       for (int i = 1; i < max; i++) {
 560         rere = WRe[k] * tem1Re[i];
 561         imim = WIm[k] * tem1Im[i];
 562         reim = WRe[k] * tem2Im[i];
 563         imre = WIm[k] * tem2Re[i];
 564
 565         temRe[radix - j] += rere + imim;
 566         temIm[radix - j] += reim - imre;
 567         temRe[j] += rere - imim;
 568         temIm[j] += reim + imre;
 569
 570         k = k + j;
 571         if (k >= radix)
 572           k = k - radix;
 573       }
 574     }
 575     for (j = 1; j < max; j++) {
 576       temRe[0] = temRe[0] + tem1Re[j];
 577       temIm[0] = temIm[0] + tem2Im[j];
 578     }
 579   }   // End of function fftPrime().
 580 }