/* thread_fft1.c  run 4, hopefully in parallel, fft1test.c */
/* see later, one file, versions  time_mp2.c  time_mp4.c   */

#include <pthread.h>
#include <string.h>	/* for strerror() */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "fft.h"

#define NTHREADS 4

#define errexit(code,str)                          \
  fprintf(stderr,"%s: %s\n",(str),strerror(code)); \
  exit(1);

/******** this is the thread code */
void *compute(void * arg)
{
  int myid=*(int *) arg;

  printf("computing %d\n",myid);
  fflush(stdout);

  /* fft1test  canned data, timing test */
  int i, j, k, m, n;
  FTYPE *a1, *a2, *b;
  int sgn = 1;	   /* default sign is 1	for FFT */
  FTYPE PI = 3.141592653589793238462643;
  FTYPE w;

  n = 65536;
  m = 32;
  w = 2.0*PI/(FTYPE)n;

  printf("fft1test: allocating space for 4*%d elements\n", n);
  b =  (FTYPE *) malloc( 2*n*sizeof(FTYPE) );  /* complex */
  a1 = (FTYPE *) malloc( n*sizeof(FTYPE) );    /* real */
  a2 = (FTYPE *) malloc( n*sizeof(FTYPE) );    /* real */
  printf("          running %d FFT and IFFT of length %d\n", m, n);

  for (i=0; i<m; i++)/* loop on data sets */
  {
    /* generate  b  test data, real, imag,  n items */
    for(j=0; j<n; j++)
    {
      *(b+2*j)   = cos(w*j)+cos(w*(i+1));
      *(b+2*j+1) = sin(w*j)+sin(w*(i+1));
    }
    if(i==0 || i==m-1)
    {
      printf("bd(%d,%d)=%g,%g,%g,%g \n", i, j, *(b+0), *(b+1), *(b+2), *(b+3));
    }
    fft(b, a1, a2, n, sgn);
    if(i==0 || i==m-1)
    {
      printf("bf(%d,%d)=%g,%g,%g,%g \n", i, j, *(b+0), *(b+1), *(b+2), *(b+3));
    }
    fft(b, a1, a2, n, -sgn);
    if(i==0 || i==m-1)
    {
      printf("bi(%d,%d)=%g,%g,%g,%g \n", i, j, *(b+0), *(b+1), *(b+2), *(b+3));
    }
  }
  return arg;
}

/******** this is the main thread's code */
int main(int argc,char *argv[])
{
  int worker;
  pthread_t threads[NTHREADS];                /* holds thread info */
  int ids[NTHREADS];                          /* holds thread args */
  int errcode;                                /* holds pthread error code */
  int *status;                                /* holds return code */
  double now, wnow;                           /* time CPU and wall */

  printf("thread_fft1.c running \n");

  now = (double)clock()/(double)CLOCKS_PER_SEC;
  wnow = (double)time(NULL);

  /* create the threads */
  for (worker=0; worker<NTHREADS; worker++) {
    ids[worker]=worker;
    if (errcode=pthread_create(&threads[worker],/* thread struct             */
		       NULL,                    /* default thread attributes */
		       compute,                 /* start routine             */
		       &ids[worker])) {         /* arg to routine            */
      errexit(errcode,"pthread_create");
    }
  }
  /* reap the threads as they exit */
  for (worker=0; worker<NTHREADS; worker++) {
    /* wait for thread to terminate */
    if (errcode=pthread_join(threads[worker],(void *) &status)) { 
      errexit(errcode,"pthread_join");
    }
    /* check thread's exit status and release its resources */
    if (*status != worker) {
      fprintf(stderr,"thread %d terminated abnormally\n",worker);
      exit(1);
    }
  }
  now = (double)clock()/(double)CLOCKS_PER_SEC - now;
  printf("CPU time is %f seconds \n", now);
  wnow = (double)time(NULL) - wnow;
  printf("wall time is %f seconds \n", wnow);
  printf("thread_fft1.c exiting \n");
  return(0);
} /* end thread_fft1.c */