Centre Blaise Pascal » Bench4GPU

SOURCE=FFT2D.c

CC=gcc

CFLAGS=-Wall -O3

LDFLAGS=-lm

CUDADIR=/opt/cuda

CUDASRC=$(CUDADIR)/src

CUDAINC=$(CUDADIR)/include

CUDALIB=$(CUDADIR)/lib64

ACML=/opt/acml

ACMLINC=$(ACML)/gfortran64_mp/include

ACMLLIB=$(ACML)/gfortran64_mp/lib

EXECUTABLE=fftw3 cufft

FORMAT=DOUBLE

#FORMAT=FLOAT

#DIRECTIVES=-D$(FORMAT) -DPRINT -DUNIT

#DIRECTIVES=-D$(FORMAT) -DUNIT -DRESULTS -DQUIET

DIRECTIVES=-DUNIT -DQUIET

all: $(EXECUTABLE)

fftw3: $(SOURCE)

	$(CC) $(CFLAGS) $(DIRECTIVES) -DFLOAT -DFFTW3 $(LDFLAGS) \

		$(SOURCE) \

		-lm -lfftw3f_threads -lfftw3f -lpthread \

		-o $(SOURCE:.c=)_SP_$@

	$(CC) $(CFLAGS) $(DIRECTIVES) -DDOUBLE -DFFTW3 $(LDFLAGS) \

		$(SOURCE) \

		-lm -lfftw3_threads -lfftw3 -lpthread \

		-o $(SOURCE:.c=)_DP_$@

cufft: $(SOURCE)

	$(CC) -I$(CUDAINC) -L$(CUDALIB) $(CFLAGS) -DFLOAT \

		-DCUFFT $(LDFLAGS) \

		$(DIRECTIVES) $(SOURCE) -lcufft -o $(SOURCE:.c=)_SP_$@

	$(CC) -I$(CUDAINC) -L$(CUDALIB) $(CFLAGS) -DDOUBLE \

		-DCUFFT $(LDFLAGS) \

		$(DIRECTIVES) $(SOURCE) -lcufft -o $(SOURCE:.c=)_DP_$@

clean: $(SOURCE)

	find . -name "$(SOURCE:.c=)_*" -exec rm {} \;

	find . -name "*~" -exec rm {} \;

check: $(EXECUTABLE)

	$(SOURCE:.c=)_SP_$(EXECUTABLE) 2 1

/* 

   Performs matrix multiply

   Thanks for help from aurel32@debian.org

*/

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

#include <sys/time.h>

#include <string.h>

#ifdef CUFFT

#include <complex.h>

#include <cufft.h>

#else

#include <fftw3.h>

#endif

#ifdef DOUBLE

#define LENGTH double

#else

#define LENGTH float

#endif

// kind of plan estimations for FFTW3

#define ESTIMATE 1

#define MEASURE 2

#define PATIENT 3

#define EXHAUSTIVE 4

// Default values

#define NTHREADS 1

#define NLOOPS 10

#define DIMENSION 1024

// Redefine FFTW function calls 

#ifdef FLOAT

#define FFTW_COMPLEX fftwf_complex

#define FFTW_PLAN fftwf_plan

#define FFTW_INIT_THREADS fftwf_init_threads

#define FFTW_PLAN_WITH_NTHREADS fftwf_plan_with_nthreads

#define FFTW_MALLOC fftwf_malloc

#define	FFTW_PLAN_DFT_2D fftwf_plan_dft_2d

#define FFTW_EXECUTE fftwf_execute

#define FFTW_DESTROY_PLAN fftwf_destroy_plan

#define FFTW_FREE fftwf_free

#elif DOUBLE

#define FFTW_COMPLEX fftw_complex

#define FFTW_PLAN fftw_plan

#define FFTW_INIT_THREADS fftw_init_threads

#define FFTW_PLAN_WITH_NTHREADS fftw_plan_with_nthreads

#define FFTW_MALLOC fftw_malloc

#define	FFTW_PLAN_DFT_2D fftw_plan_dft_2d

#define FFTW_EXECUTE fftw_execute

#define FFTW_DESTROY_PLAN fftw_destroy_plan

#define FFTW_FREE fftw_free

#endif

int main(int argc,char **argv)

{

  int i;

  int nloops=NLOOPS,nthreads=NTHREADS,size;

  int type_plan=0;

  double dplan=0.,dcompute=0.;

  struct timeval tv1,tv2;

  struct timezone tz;

  if ((argc==1)||

      (strcmp(argv[1],"-h")==0)||

      (strcmp(argv[1],"--help")==0))

    {

      printf("\n\tCalculate fake FFTW images\n\n"

	     "\tParameters to give :\n\n"

	     "\t1> Plan ESTIMATE MEASURE PATIENT EXHAUSTIVE\n"

	     "\t2> Size of Image (2^n)\n"

	     "\t3> Number of Loops (integer)\n"

	     "\t4> Number of Threads (integer)\n\n");

      return(0);

    }

  // initialization of FFTW threads

  if (atoi(argv[4])<1) {

    nthreads=1;

  }

  else {

    nthreads=atoi(argv[4]);

  }

  if (atoi(argv[3])<1) {

    nloops=1;

  }

  else {

    nloops=atoi(argv[3]);

  }

  if (atoi(argv[2])<1) {

    size=2;

  }

  else {

    size=atoi(argv[2]);

  }

  printf("%i %i\n",size,nthreads);

#ifdef CUFFT

  cufftHandle plan;

  cufftComplex *in, *devin;

  size_t arraySize = sizeof(cufftComplex)*size*size;

  cudaMallocHost((void**) &in, arraySize);

  cudaMalloc((void**) &devin, arraySize);

  // initialisation of arrays

  for (i = 0; i < size*size; i++) {

    in[i].x=1.;

    in[i].y=0.;

  }

  in[0].x=1.;

  in[0].y=0.;

  // First timer to start plan

  gettimeofday(&tv1, &tz);

  // Plan & copy to device

  cufftPlan2d(&plan, size, size, CUFFT_C2C);

  cudaMemcpy(devin, in, arraySize, cudaMemcpyHostToDevice);

  // Second timer to end plan      

  gettimeofday(&tv2, &tz);	

  dplan=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +		\

		     (tv2.tv_usec-tv1.tv_usec))/1000000.;  

  // First timer to start process

  gettimeofday(&tv1, &tz);

  // Process

  for (i=0;i<nloops;i++) {

    cufftExecC2C(plan, devin, devin, CUFFT_FORWARD);

  }

  cudaMemcpy(in, devin, arraySize, cudaMemcpyDeviceToHost);

  printf("%i=%f,%f\n",0,in[0].x,in[0].y);

  // Second timer to end process

  gettimeofday(&tv2, &tz);

  cufftDestroy(plan);

  cudaFreeHost(in);

  cudaFree(devin);

#elif FFTW3      

  FFTW_COMPLEX *in;

  FFTW_PLAN p;

  FFTW_INIT_THREADS();

  FFTW_PLAN_WITH_NTHREADS(nthreads);

  in = (FFTW_COMPLEX*) FFTW_MALLOC(sizeof(FFTW_COMPLEX)*size*size);

  // First timer to start plan

  gettimeofday(&tv1, &tz);

  if (strcmp(argv[1],"MEASURE")==0) {

    p = FFTW_PLAN_DFT_2D(size,size,in,in,FFTW_FORWARD,FFTW_MEASURE);

    type_plan=MEASURE;

  }

  if (strcmp(argv[1],"PATIENT")==0) {

    p = FFTW_PLAN_DFT_2D(size,size,in,in,FFTW_FORWARD,FFTW_PATIENT);

    type_plan=PATIENT;

  }

  if (strcmp(argv[1],"EXHAUSTIVE")==0) {

    p = FFTW_PLAN_DFT_2D(size,size,in,in,FFTW_FORWARD,FFTW_EXHAUSTIVE);

    type_plan=EXHAUSTIVE;

  }

  if ((type_plan==0)||(strcmp(argv[1],"ESTIMATE")==0)) {

    p = FFTW_PLAN_DFT_2D(size,size,in,in,FFTW_FORWARD,FFTW_ESTIMATE);

    type_plan=ESTIMATE;

  }

  // Second timer to end plan      

  gettimeofday(&tv2, &tz);	

  dplan=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +		\

		 (tv2.tv_usec-tv1.tv_usec))/1000000.;  

  // initialisation of arrays

  for (i = 0; i < size*size; i++) {

    (in[i])[0]=0.;

    (in[i])[1]=0.;

  }

  (in[0])[0]=1.;

  (in[0])[1]=0.;

  // First timer to start process

  gettimeofday(&tv1, &tz);

  for (i=0;i<nloops;i++) {

    FFTW_EXECUTE(p);

  }

  // Second timer to end process

  gettimeofday(&tv2, &tz);

  FFTW_DESTROY_PLAN(p);

  FFTW_FREE(in); 

#endif

  dcompute=(double)((tv2.tv_sec-tv1.tv_sec) * 1000000L +		\

		    (tv2.tv_usec-tv1.tv_usec))/1000000./(double)nloops;  

  printf("%s,%i,%i,%i,%i,%lf,%lf\n",argv[1],

	 size,size,nthreads,nloops,dplan,dcompute);

  /* if ((file=fopen(argv[7],"a"))==NULL) */

  /* 	{ */

  /* 	  printf("Impossible to append result in %s file\n", */

  /* 		 argv[7]); */

  /* 	  return(0); */

  /* 	} */

  /* fprintf(file,"%s,%i,%i,%i,%i,%lf,%lf\n", */

  /* 	  argv[1],size,size,nthreads,nloops,dplan,dcompute); */

  /* fclose(file); */

  return 0;

}