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