#include "spfft.h"

/* #define expertcall*/
 
void output(arr,pos,ctrl)
struct FFTCTL *ctrl;
PLCPLX* arr;
int pos;
{
   int tadr,madr,padr;
   plural unsigned itmp1, itmpM1;
   plural FLTYP tmp0,tmp1;
  
   tadr = pos;
   mp_adr(tadr, madr, padr, ctrl->mb);
  
   tmp0 = arr[madr][0];
   tmp1 = arr[madr][1];
   printf("total address: %6d arr. contains: %8.2f %8.2f\n",
         tadr,proc[padr].tmp0 ,proc[padr].tmp1);
}

 
/* initialize data for whole array */
void init(arr, ctrl)
PLCPLX *arr;
struct FFTCTL *ctrl;
{
   int *mb;
   plural int procadr, tadr;
   int i0,i1,i2, madr;
 
   mb = ctrl->mb;
   madr = 0;
   procadr = iyproc<<(mb[1]+mb[0]+lnxproc) | (ixproc<<mb[0]);
 
   for (i2=0; i2< (1<<(mb[2]+1)); i2++) {
      for (i1=0; i1< (1<<mb[1]); i1++) {
         for (i0=0; i0< (1<<mb[0]); i0++) {
            tadr = (i2<<(mb[1]+mb[0]+lnproc)) | (i1<<(lnxproc+mb[0]))
                   | mb[0] | procadr;
            arr[madr][0] = p_sqrt((plural FLTYP) tadr);
            arr[madr][1] = 0;
            madr++;
         }
      }
   }
}
 

main()
{
   struct FFTCTL *ctrl;
   PLCPLX *arr;
   plural unsigned itmp1,itmpM1;
   plural FLTYP tmp0,tmp1;
   unsigned madr,asiz,tadr,padr;
   int trlng, memorg, nb;
   int igr,ngr,mb[3],mbts,agrlst[64];
   int asizx,asizy,trlngx,trlngy,numdim;

   printf("Arraysize x-dir, arraysize y-dir: ");
   scanf("%d%d", &asizx, &asizy);
   asiz = asizx*asizy/nproc;
   P_MALLOC(arr,asiz,PLCPLX);
   numdim = 2;
   printf("Give transform length x, transform length y: ");
   scanf("%d%d", &trlngx, &trlngy);
   ctrl = fft_init(1,numdim,asizx,trlngx,asizy,trlngy,0,0);
   
   /* initialize data */
   init(arr, ctrl);

   printf("here are some input\n");
   output(arr,0,ctrl);
   output(arr,1,ctrl);
   output(arr,457,ctrl);
   output(arr,16383,ctrl);
   output(arr,32000,ctrl);
   output(arr,asiz*nproc-1,ctrl);
   
   /* Check code by running forw. fft followed by inverse fft
      ordered -> ordered -> ordered */
   /* do forwards fft */
   sp_fft(arr,1,ctrl); 

   printf("here are some output, after sorted fft\n");
   output(arr,0,ctrl);
   output(arr,1,ctrl);
   output(arr,457,ctrl);
   output(arr,16383,ctrl);
   output(arr,32000,ctrl);
   output(arr,asiz*nproc-1,ctrl);

   /* do inverse fft */
   sp_fft(arr,-1,ctrl);
   

   printf("here are some output (after inv.fft+scale) \n"); 
   output(arr,0,ctrl);
   output(arr,1,ctrl);
   output(arr,457,ctrl);
   output(arr,16383,ctrl);
   output(arr,32000,ctrl);
   output(arr,asiz*nproc-1,ctrl);

   {
   printf("Timing:\n");
   printf(
      "Flop rate are based on the cost of a similar complex transform\n");
   {
      int nops,lg2n,numiter,numit,i,ib;
      unsigned ibmsk;
      struct timeval beg,end;
      double stime, etime, eltime;
      float mflop;
 
      /* count no. of active bits */
      for (igr=0,ibmsk=1,lg2n=0;igr<32;igr++,ibmsk<<=1)
         if ((ctrl->actbits) & ibmsk) lg2n++;
      nops = 5*asiz*nproc*lg2n;
      numiter = 2*(300000000)/nops + 1;
 
      /* initialize data */
      for (madr=0;madr<asiz;madr++) {
         arr[madr][0] = 0;
         arr[madr][1] = 0;
      }
 
      numit = numiter/10 + 1;
      gettimeofday(&beg,0);
      {
         fft_free(ctrl);
#ifndef expertcall
         ctrl = fft_init(1,numdim,asizx,trlngx,asizy,trlngy,0,0);
         /* first call to fft_srt and sp_fft does some initializations */
         sp_fft(arr,1,ctrl);
         sp_fft(arr,2,ctrl);
#else
         ctrl = fft_ainit(1,mb,ngr,agrlst,NULL,NULL);
         /* first call to fft_srt and sp_fft does some initializations */
         sp_fft(arr,1,ctrl);
         sp_fft(arr,2,ctrl);
#endif
      }
      gettimeofday(&end,0);
      stime = beg.tv_sec + 0.000001*beg.tv_usec;
      etime = end.tv_sec + 0.000001*end.tv_usec;
      eltime = (etime - stime);
      mflop = ((float) nops*numit)/(1000000*eltime);
      /* printf("Initialization speed:  %6.2f'megainit'\n",mflop); */
 
      gettimeofday(&beg,0);
      for (i=0;i<numiter;i++) {
         sp_fft(arr,2,ctrl);
      }
      gettimeofday(&end,0);
      stime = beg.tv_sec + 0.000001*beg.tv_usec;
      etime = end.tv_sec + 0.000001*end.tv_usec;
      eltime = (etime - stime);
      mflop = ((float) nops*numiter)/(1000000*eltime);
      printf("Forward transforms\n");
      printf("   Speed without unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(numiter));

      gettimeofday(&beg,0);
      for (i=0;i<numiter;i++) {
         sp_fft(arr,-2,ctrl);
      }
      gettimeofday(&end,0);
      stime = beg.tv_sec + 0.000001*beg.tv_usec;
      etime = end.tv_sec + 0.000001*end.tv_usec;
      eltime = (etime - stime);
      mflop = ((float) nops*numiter)/(1000000*eltime);
      printf("Inverse transforms\n");
      printf("   Speed without unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(numiter));
 
      gettimeofday(&beg,0);
      for (i=0;i<numiter;i++) {
         sp_fft(arr,1,ctrl);
      }
      gettimeofday(&end,0);
      stime = beg.tv_sec + 0.000001*beg.tv_usec;
      etime = end.tv_sec + 0.000001*end.tv_usec;
      eltime = (etime - stime);
      mflop = ((float) nops*numiter)/(1000000*eltime);
      printf("Forward transforms\n");
      printf("   Speed with unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(numiter));

      gettimeofday(&beg,0);
      for (i=0;i<numiter;i++) {
         sp_fft(arr,-1,ctrl);
      }
      gettimeofday(&end,0);
      stime = beg.tv_sec + 0.000001*beg.tv_usec;
      etime = end.tv_sec + 0.000001*end.tv_usec;
      eltime = (etime - stime);
      mflop = ((float) nops*numiter)/(1000000*eltime);
      printf("Inverse transforms\n");
      printf("   Speed with unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(numiter));
   }
   }
}