example3.cpp

基于GPU进行快速科学计算

/******************************************************************************\

  Copyright 2005 The University of North Carolina at Chapel Hill.
  All Rights Reserved.

  Permission to use, copy, modify and distribute this software and its
  documentation for educational, research and non-profit purposes, without
  fee, and without a written agreement is hereby granted, provided that the
  above copyright notice and the following three paragraphs appear in all
  copies. Any use in a commercial organization requires a separate license.

  IN NO EVENT SHALL THE UNIVERSITY OF NORTH CAROLINA AT CHAPEL HILL BE LIABLE
  TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
  DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND
  ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF NORTH CAROLINA HAVE BEEN
  ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.


  Permission to use, copy, modify and distribute this software and its
  documentation for educational, research and non-profit purposes, without
  fee, and without a written agreement is hereby granted, provided that the
  above copyright notice and the following three paragraphs appear in all
  copies.

  THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES,
  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
  FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN
  "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO
  PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.


   ---------------------------------
  |Please send all BUG REPORTS to:  |
  |                                 |
  |   [e-mail]               |
  |                                 |
   ---------------------------------


  The authors may be contacted via:

  US Mail:         N. Govindaraju 
                       Department of Computer Science
                       Sitterson Hall, CB #3175
                       University of North Carolina
                       Chapel Hill, NC 27599-3175

\*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "stopwatch.hpp"
#include <GPUFFTW.h>
#include <math.h>
#define ITERATIONS 1


#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

//typedef double FFT_DATA_TYPE ;
#define FFT_DATA_TYPE float

//#define REAL_FFT 1


int main(int argc, char *argv[]){
	Stopwatch gputimer, cputimer;	//timers 
	FFT_ERROR error;	//sorting error code
	float *array=NULL;	//array to be passed to the GPU for sorting
	float *test=NULL;	//array to be used for CPU sorting
	GPUFFTW *gpufft;	//GPU sorting object

	//Precision to use
	printf("\nTesting in 32-bit precision\n");
		
	//Create object for sorting elements in 32-bit specified precision
	gpufft = new GPUFFTW();			
		
	//Find the maximum array size supported by the video card
	int S = gpufft->maxSize();
	array=NULL;
	//Find out how big a set of arrays the system RAM can support
	while(array==NULL){
		array=new float[4*S];
		S=S/2;
	}
	S=S*2;
	delete[] array;
	float  *imag = new float[S];
	float *real = new float[S];

	FFT_DATA_TYPE *cpureal = new FFT_DATA_TYPE[S];
	FFT_DATA_TYPE *tmpreal = new FFT_DATA_TYPE[S];

	FFT_DATA_TYPE *cpuimag = new FFT_DATA_TYPE[S];
	FFT_DATA_TYPE *tmpimag = new FFT_DATA_TYPE[S];

	FFT_DATA_TYPE *imagdata = new FFT_DATA_TYPE[S];
	FFT_DATA_TYPE *realdata = new FFT_DATA_TYPE[S];


	
	printf("\n\nYour Video Card can sort arrays of sizes upto: %d",gpufft->maxSize());
	printf("\nMaximum input size to be tested on your system (determined by amount of system RAM): %d",S);
	printf("\nPRECISION: %d-bit\n",32);
	
	int logsize=9;
	//Try different input array sizes
	for(int size=(1<<logsize);size<=(1<<22);size*=2, logsize++)
	{
		printf("Input Sequence Size: %d\n",size);
		gputimer.Reset();
		cputimer.Reset();		
		int t=0;
		for(int i=0;i<=ITERATIONS;i++){
			//generate records with random numbers
			for(int j=0;j<size;j++){
				real[j] = (((float)rand())/RAND_MAX)-0.5;
				imag[j] = (((float)rand())/RAND_MAX)-0.5;
			}

			//Make a copy of input
			
			int d;
			
			for(d=0; d<size;d++){
				cpureal[d] = real[d];
				tmpreal[d] = real[d];
				cpuimag[d] = imag[d];
				tmpimag[d] = imag[d];
				realdata[d]= real[d];
				imagdata[d] = imag[d];
			}
			/*memcpy(cpureal,real,size*sizeof(float));
			memcpy(cpuimag,imag,size*sizeof(float));
			memcpy(tmpreal,real,size*sizeof(float));
			memcpy(tmpimag,imag,size*sizeof(float));
			*/
			if(i) gputimer.Start();
#ifndef REAL_FFT
			error=gpufft->cgpufft1d(real,imag,size,1.0);
#else
			error=gpufft->scgpufft1d(real,size,1.0);
#endif
			if(i) gputimer.Stop();		
			

			
			if(error!=FFT_SUCCESS){
				char *error_string = gpufft->FFTErrorString(error);
				printf("\n%s\n",error_string);
				delete[] error_string;
				delete[] real;
				delete[] imag;
				delete(gpufft);
				return(1);
			}

#if 0		//use it to test the precision and correctness of the algorithm	
			//fft code on cpu
			int count=0;
			if(i) cputimer.Start();
			for(int t=1; t<=logsize;t++){
				for(int fftj = 0; fftj< (1<< (logsize-t)); fftj++){
					for(int fftk=0; fftk< (1<< (t-1)); fftk++){
						int index1 = fftj * (1<<t)+fftk;
						int index11 = fftj * (1<<(t-1))+fftk;
						int index12 = fftj * (1<<(t-1))+ size/2 + fftk;

						int index2 = index1 + (1<<(t-1));
						double angle = -2.0 * (M_PI) * fftk / (1<<t);
						
						tmpreal[index1] = cpureal[index11] + cos(angle) * cpureal[index12] - sin(angle) *cpuimag[index12];
						tmpimag[index1] = cpuimag[index11] + cos(angle) * cpuimag[index12] + sin(angle) *cpureal[index12];

						tmpreal[index2] = cpureal[index11] -(cos(angle) * cpureal[index12] - sin(angle) *cpuimag[index12]);
						tmpimag[index2] = cpuimag[index11] -(cos(angle) * cpuimag[index12] + sin(angle) *cpureal[index12]);
						
						
					}
				}
				
				memcpy(cpureal,tmpreal,size*sizeof(FFT_DATA_TYPE));
				memcpy(cpuimag,tmpimag,size*sizeof(FFT_DATA_TYPE));

			}
			if(i) cputimer.Stop();
#endif
#ifndef REAL_FFT			
			error=gpufft->cgpufft1d(real, imag,size,-1);
#else
			error=gpufft->scgpufft1d(real, size,-1);
#endif
			

			double diff=0.0;
			double magnitude=0.0;

			for(d=0; d<size;d++){
#ifndef REAL_FFT
				diff+=fabs(imagdata[d] - imag[d])*fabs(imagdata[d] - imag[d]);
#endif
				diff+=fabs(realdata[d] - real[d])*fabs(realdata[d] - real[d]);
#ifndef REAL_FFT
				magnitude+= real[d]*real[d]+imag[d]*imag[d];
#else
				magnitude+= real[d]*real[d];
#endif

			//	if(   fabsf(imagdata[d] - imag[d])> 0.0015 ||1) printf("ERROR: IMAG: %f %f : %d\n", imagdata[d], imag[d], d);
			//	if(   fabsf(cpureal[d] -  real[d])> 0.0015 ||1) printf("ERROR: REAL: %f %f : %d\n", cpureal[d], real[d], d);
			}
			double precision = sqrt(diff)/sqrt(magnitude);
			printf("precision (L2):%4.10f \n", precision);
			
			
		}
		printf(" GPUFFTW Time: ");			
		printf("%.3f sec",(gputimer.GetTime())/ITERATIONS);
		printf(" Qsort() Time: ");
		printf("%.3f sec",(cputimer.GetTime())/ITERATIONS);
		printf("  DIFFERENCE: %.3f sec\n",(float)cputimer.GetTime()-gputimer.GetTime());
		

	}
	
	return(0);
}