#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]); 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;

#ifndef expertcall 
   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(0,numdim,asizx,trlngx,asizy,trlngy,0,0);
#else
   printf("Memory org.: mb[0],mb[1],mb[2] : ");
   scanf("%d%d%d",&mb[0],&mb[1],&mb[2]);
   mbts = mb[0]+mb[1]+mb[2];
   asiz = 1<<mbts;
   P_MALLOC(arr,asiz,PLCPLX);
   printf("No. of active groups: ");
   scanf("%d",&ngr);
   for(igr=0;igr<ngr;igr++) {
      printf("Group %d. Enter startbit and length (in bits) : ",igr);
      scanf("%d%d",&agrlst[2*igr],&agrlst[2*igr+1]);
   }
   ctrl = fft_ainit(0,mb,ngr,agrlst,NULL,NULL);
#endif
   
   /* 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 -> scrambled -> ordered */
   /* do forwards fft */
   sp_fft(arr,1,ctrl); 

   printf("Here are some sorted fourier coeficients\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);

   /* initialize data */
   init(arr, ctrl);

   /* Check code by running forw. fft + sorting + scaling 4 times */
   /* do forwards fft */
   sp_fft(arr,1,ctrl);
   sp_fft(arr,1,ctrl);
   /* scale the data */
   {
   FLTYP scalfac;
   scalfac = ctrl->scale;
   for (madr=0;madr<asiz;madr++) {
      arr[madr][0] *= scalfac;
      arr[madr][1] *= scalfac;
   }
   }
   /* do forwards fft */
   sp_fft(arr,1,ctrl);
   sp_fft(arr,1,ctrl);
   /* scale the data */
   {
   FLTYP scalfac;
   scalfac = ctrl->scale;
   for (madr=0;madr<asiz;madr++) {
      arr[madr][0] *= scalfac;
      arr[madr][1] *= scalfac;
   }
   }
   printf("here are some output (after 4 times sorted forw. 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);

   /* SOME TIMINGS: */
   {
   printf("Timing:\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 = (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);
      for (i=0;i<numit;i++) { 
         fft_free(ctrl);
#ifndef expertcall
         ctrl = fft_init(0,numdim,asizx,trlngx,asizy,trlngy,0,0);
         /* first call to fft_srt and sp_fft does some initializations */
         fft_srt(arr,0,ctrl);
         sp_fft(arr,0,ctrl);
#else
         ctrl = fft_ainit(0,mb,ngr,agrlst,NULL,NULL);
         /* first call to fft_srt and sp_fft does some initializations */
         fft_srt(arr,0,ctrl);
         sp_fft(arr,0,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);
         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 = 2*((float) nops*numiter)/(1000000*eltime);
      printf("Speed without unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(2*numiter));

      gettimeofday(&beg,0);
      for (i=0;i<numiter;i++) { 
         sp_fft(arr,1,ctrl);
         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 = 2*((float) nops*numiter)/(1000000*eltime);
      printf("Speed with unscrambling: %6.2fMflop\n",mflop);
      printf("   Elapsed time: %f millisec. \n",eltime*1000/(2*numiter));
   }
   } 
}