/* ---------------------------------------------------------------
 * Serial - Two-Dimensional Fast Fourier Transform - C Version
 * FILE: ser_2dfft.c
 * OTHER FILES: ser_fft.c ser_fft.h
 * DESCRIPTION: Algorithm: Perform a 1D FFT on each row, transpose, 
 * then perform a 1D FFT on each row again.  
 *
 * Input is an 512x512 complex matrix. Output is an 512x512
 * complex matrix that overwrites the input matrix.  The input matrix is
 * initialized with a point source.
 *
 * A straightforward unsophisticated 1D FFT kernel is used.  It is
 * sufficient to convey the general idea, but be aware that there are
 * better 1D FFTs available on many systems.
 *
 * AUTHOR: George Gusciora
 * LAST REVISED:  
 * --------------------------------------------------------------- */

#include "ser_fft.h"

#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#define IMAGE_SIZE              512

static mycomplex a[MAXN][MAXN];     /* input matrix */
static mycomplex w_common[MAXN/2];  /* twiddle factors */

main(argc,argv)
int argc;
char **argv;
{
  int n;           /* FFT size */
  int logn;        /* log base 2 of n */
  int errors,sign; /* used for error checking */
  int nx;          /* used to compute logn */
  int flops;       /* total number of floating point ops */
  float mflops;    /* Mflops/s */
  double fsecs;    /* time spent doing 1D FFT's */
  double tsecs;    /* time spent doing the transpose */
  double secs;     /* total time of the 2D FFT */
  int i,j;         /* index variables */

  struct timeval start, finish; /* gettimeofday structures */
  void print_cmat(), print_cvec(); /* forward definitions */


  n = IMAGE_SIZE;
  /* compute logn and ensure that n is a power of two */
  nx = n;
  logn = 0;
  while(( nx >>= 1) > 0) 
    logn++; 
  nx = 1;
  for (i=0; i<logn; i++)
    nx = nx*2;
  if (nx != n) {
    (void)fprintf(stderr, "%s: fft size must be a power of 2\n",argv[0]);
    exit(0);
  }
  
  /* initialize the input matrix with a centered point source */
  for (i=0; i<n; i++)
    for (j=0; j<n; j++) 
      a[i][j].r = a[i][j].i = 0.0;
  a[n/2][n/2].r =  a[n/2][n/2].i = (float)n; 

  /* precompute the complex constants (twiddle factors) for the 1D FFTs */
  for (i=0;i<n/2;i++) {
    w_common[i].r = (float) cos((double)((2.0*PI*i)/(float)n));
    w_common[i].i = (float) -sin((double)((2.0*PI*i)/(float)n));
  }
  
  /* 
   * now do the actual 2D FFT
   */

  /* first do a set of row FFTs */
  (void)gettimeofday(&start, 0);
  for (i=0; i<n; i++)
    (void)fft(&a[i][0], w_common, n, logn);
  (void)gettimeofday(&finish, 0);
  fsecs = ((((finish.tv_sec - start.tv_sec) * 1000000.0) +
	    finish.tv_usec) - start.tv_usec) / 1000000.0;

  /* then transpose the matrix */
  (void)gettimeofday(&start, 0);
  for (i=0; i<n; i++)
    for (j=i+1; j<n; j++) 
      SWAP(a[i][j], a[j][i]); 
  (void)gettimeofday(&finish, 0);
  tsecs = ((((finish.tv_sec - start.tv_sec) * 1000000.0) +
	    finish.tv_usec) - start.tv_usec) / 1000000.0;
      
  /* then do another set of row FFTs */
  (void)gettimeofday(&start, 0);
  for (i=0; i<n; i++)
    (void)fft(&a[i][0], w_common, n, logn);
  (void)gettimeofday(&finish, 0);
  fsecs += ((((finish.tv_sec - start.tv_sec) * 1000000.0) +
	    finish.tv_usec) - start.tv_usec) / 1000000.0;

  /* transpose the matrix back */
  (void)gettimeofday(&start, 0);
  for (i=0; i<n; i++)
    for (j=i+1; j<n; j++) 
      SWAP(a[i][j], a[j][i]); 
  (void)gettimeofday(&finish, 0);
  tsecs = ((((finish.tv_sec - start.tv_sec) * 1000000.0) +
	    finish.tv_usec) - start.tv_usec) / 1000000.0;
      
/*  print_cmat(a, n, n);*/

  /* check the answers for an alternating sequence of (+-n,+-n) */
  errors = 0;
  for (i=0; i<n; i++) {
    if (((i+1)/2)*2 == i) 
      sign = 1;
    else
      sign = -1;
    for (j=0; j<n; j++) {
      if (a[i][j].r > n*sign+EPSILON ||
	  a[i][j].r < n*sign-EPSILON ||
	  a[i][j].i > n*sign+EPSILON ||
	  a[i][j].i < n*sign-EPSILON) {
	printf("%d,%d\n",i,j);
	errors++;
      }
      sign *= -1;
    }
  }
  if (errors) { 
    printf("%d errors!!!!!\n", errors);
    exit(0);
  }

  /* summarize the 2d FFT performance */
  printf("%d x %d 2D FFT\n", n, n);
  secs = fsecs + tsecs;
  flops = (n*n*logn)*10;
  mflops = ((float)flops/1000000.0);
  mflops = mflops/(float)secs;
  printf("\n");
  printf("1D FFTs   : %10.6f secs (%2d%%)\n", fsecs, (int)(fsecs/secs*100));
  printf("transpose : %10.6f secs (%2d%%)\n", tsecs, (int)(tsecs/secs*100));
  printf("total     : %10.6f secs\n", secs);
  printf("\n");
  printf("%10.6f Mflop/s\n", mflops);

  exit(0);
}

/*
 * routines to print complex matrices and vectors
 */
void print_cvec(a, n)
mycomplex *a;
int n;
{
  int i;

  for (i=0; i<n; i++) {
    if (i%4 == 0 && i)  
      printf("\n"); 
    printf("(%.3f,%.3f) ", a[i].r, a[i].i);   
  }
  printf("\n");
}

void print_cmat(a, m, n)
mycomplex a[][MAXN];
int m,n;
{
  int i,j;

  for (i=0; i<m; i++) {
    for (j=0; j<n; j++)
      printf("(%6.2f,%6.2f) ",a[i][j].r, a[i][j].i);
    printf("\n");
  }
}