fft.cpp

fft并行算法实现

// Fft.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mpi.h"

#define MAX_N 50
#define PI 3.141592654
#define EPS 10E-8
#define V_TAG 99
#define P_TAG 100
#define Q_TAG 101
#define R_TAG 102
#define S_TAG 103
#define S_TAG2 104

//typedef enum {false,true}bool;
typedef struct
{
	double r;
	double i;
}complex_t;

complex_t p[MAX_N],q[MAX_N],s[2*MAX_N],r[2*MAX_N];
complex_t w[2*MAX_N];
int variableNum;
double transTime=0,totalTime=0,beginTime;
MPI_Status status;

/**
   函数，复数加
**/
void comp_add(complex_t* result,const complex_t *c1,const complex_t *c2)
{
	result->i=c1->i+c2->i;
	result->r=c1->r+c2->r;
}

/**
   函数，复数乘
**/
void comp_multiply(complex_t* result,const complex_t *c1,const complex_t *c2)
{
	result->i=c1->r*c2->i+c1->i*c2->r;
	result->r=c1->r*c2->r-c1->i*c2->i;
}

/**
   函数，数组重新排列
**/
void shuffle(complex_t * f,int beginpos,int endpos)
{
	int i,j;
	complex_t temp[2*MAX_N];

	for(i=beginpos;i<=endpos;i++)
	{
		temp[i]=f[i];
	}

	j=beginpos;
	for(i=beginpos;i<=endpos;i+=2)
	{
		f[j]=temp[i];
		j++;
	}
	for(i=beginpos+1;i<=endpos;i+=2)
	{
		f[j]=temp[i];
		j++;
	}
}


/**
   函数：对复数序列进行FFT或者IFFT
**/
void evaluate(complex_t * f,int beginpos,int endpos,
			  const complex_t * x, complex_t * y,int leftpos,int rightpos,int totallength)
{
	int i;
	if((beginpos>endpos)||(leftpos>rightpos))
		exit(-1);
	else if(beginpos == endpos)
	{
		for(i=leftpos;i<=rightpos;i++)
			y[i]=f[beginpos];
	}
	else if(beginpos+1==endpos)
	{
		for(i=leftpos;i<=rightpos;i++)
		{
			complex_t temp;
			comp_multiply(&temp,&f[endpos],&x[i]);
			comp_add(&y[i],&f[beginpos],&temp);
		}
	}
	else
	{
		complex_t tempX[2*MAX_N];
		complex_t tempY1[2*MAX_N],tempY2[2*MAX_N];
		int midpos=(beginpos+endpos)/2;
		shuffle(f,beginpos,endpos);

		for(i=leftpos;i<=rightpos;i++)
			comp_multiply(&tempX[i],&x[i],&x[i]);
		evaluate(f,beginpos,midpos,tempX,tempY1,leftpos,rightpos,totallength);
		evaluate(f,midpos+1,endpos,tempX,tempY2,leftpos,rightpos,totallength);

		for(i=leftpos;i<=rightpos;i++)
		{
			complex_t temp;
			comp_multiply(&temp,&x[i],&tempY2[i]);
			comp_add(&y[i],&tempY1[i],&temp);
		}
	}
}

/**
   函数：打印元素实部
**/
void print(const complex_t *f,int flength)
{
	bool isprint=false;
	int i;

	if((f[0].r>0 && f[0].r>EPS)||(f[0].r<0 && f[0].r<-EPS))
	{
		printf("%f",f[0].r);
		isprint=true;
	}
	for(i=1;i<flength;i++)
	{
		if(f[i].r>EPS)
		{
			if(isprint)
				printf("+");
			else
				isprint=true;
			printf("%ft^%d",f[i].r,i);
		}
		else if(f[i].r<-EPS)
		{
			if(isprint)
				printf("-");
			else
				isprint=true;
			printf("%ft^%d",-f[i].r,i);
		}
	}
	if(isprint==false)
		printf("0");
	printf("\n");
}

/**
   函数：完整打印数组元素
**/
void myprint(const complex_t* f,int flength)
{
	int i;
	for(i=0;i<flength;i++)
	{
		printf("%f+%fi, ",f[i].r,f[i].i);
	}
	printf(" \n");
}

/**
   函数：累计时间
**/
void addTransTime(double t)
{
	transTime+=t;
}

/**
   函数：从文件中读取数据
**/
bool readFromFile()
{
	int i;
	FILE* fin=fopen("data.txt","r");
	if(fin==NULL)
	{
		printf("cannot open data file-data.txt\n");
		return false;
	}
	fscanf(fin,"%d\n",&variableNum);
	if((variableNum<1)||(variableNum>MAX_N))
	{
		printf("variablenum out of range\n");
		return false;
	}

	for(i=0;i<variableNum;i++)
	{
		fscanf(fin,"%lf",&p[i].r);
		p[i].i=0;
	}
	for(i=0;i<variableNum;i++)
	{
		fscanf(fin,"%lf",&q[i].r);
		q[i].i=0;
	}

	fclose(fin);

	printf("read from data file\n");
	printf("p(t)=");
	print(p,variableNum);
	printf("q(t)=");
	print(q,variableNum);

	return true;
}

/**
   函数：发送原始数据
**/
void sendOrigData(int size)
{
	int i;
	for(i=1;i<size;i++)
	{
		MPI_Send(&variableNum,1,MPI_INT,i,
			V_TAG,MPI_COMM_WORLD);
		MPI_Send(p,variableNum*2,MPI_DOUBLE,i,
			P_TAG,MPI_COMM_WORLD);
		MPI_Send(q,variableNum*2,MPI_DOUBLE,i,
			Q_TAG,MPI_COMM_WORLD);
	}
}

/**
   函数：接收原始数据
**/
void recvOrigData()
{
	MPI_Recv(&variableNum,1,MPI_INT,0,V_TAG,MPI_COMM_WORLD,&status);
	MPI_Recv(p,variableNum*2,MPI_DOUBLE,0,P_TAG,MPI_COMM_WORLD,&status);
	MPI_Recv(q,variableNum*2,MPI_DOUBLE,0,Q_TAG,MPI_COMM_WORLD,&status);
}

int main(int argc,char * argv[])
{
	int rank,size,i;
	int wlength;
	int averagelength,morelength,startpos,stoppos;

	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
	MPI_Comm_size(MPI_COMM_WORLD,&size);

	//main computationcode here
	if(rank==0)
	{
		if(!readFromFile())
			exit(-1);
		if(size>2*variableNum)
		{
			printf("too many processors, reduce your -np value\n");
			MPI_Abort(MPI_COMM_WORLD,1);
		}
		beginTime=MPI_Wtime();
		sendOrigData(size);
		addTransTime(MPI_Wtime()-beginTime);
	}
	else
	{
		recvOrigData();
	}

	wlength=2*variableNum;
	for(i=0;i<wlength;i++)
	{
		w[i].r=cos(i*2*PI/wlength);
		w[i].i=sin(i*2*PI/wlength);
	}

	averagelength=wlength/size;
	morelength=wlength%size;
	startpos=morelength+rank*averagelength;
	stoppos=startpos+averagelength-1;
	if(rank==0)
	{
		startpos=0;
		stoppos=morelength+averagelength-1;
	}
	evaluate(p,0,variableNum-1,w,s,startpos,stoppos,wlength);
	evaluate(q,0,variableNum-1,w,r,startpos,stoppos,wlength);

	for(i=startpos;i<=stoppos;i++)
	{
		complex_t temp;
		comp_multiply(&temp,&s[i],&r[i]);
		s[i]=temp;
		s[i].r/=wlength*1.0;
		s[i].i/=wlength*1.0;
	}

	if(rank>0)
	{
		MPI_Send(s+startpos,averagelength*2,MPI_DOUBLE,0,S_TAG,MPI_COMM_WORLD);
		MPI_Recv(s,wlength*2,MPI_DOUBLE,0,S_TAG2,MPI_COMM_WORLD,&status);
	}
	else
	{
		double tempTime=MPI_Wtime();
		for(i=1;i<size;i++)
		{
			MPI_Recv(s+morelength+i*averagelength,averagelength*2,MPI_DOUBLE,i,S_TAG,MPI_COMM_WORLD,&status);
		}
		for(i=1;i<size;i++)
		{
			MPI_Send(s,wlength*2,MPI_DOUBLE,i,S_TAG2,MPI_COMM_WORLD);
		}
		addTransTime(MPI_Wtime()-tempTime);
	}

	for(i=1;i<wlength/2;i++)
	{
		complex_t temp;
		temp=w[i];
		w[i]=w[wlength-i];
		w[wlength-i]=temp;
	}
	evaluate(s,0,wlength-1,w,r,startpos,stoppos,wlength);

	if(rank>0)
	{
		MPI_Send(r+startpos,averagelength*2,MPI_DOUBLE,0,R_TAG,MPI_COMM_WORLD);
	}
	else
	{
		double temptime=MPI_Wtime();
		for(i=1;i<size;i++)
		{
			MPI_Recv((r+morelength+i*averagelength),averagelength*2,MPI_DOUBLE,i,R_TAG,
				MPI_COMM_WORLD,&status);
		}
		totalTime=MPI_Wtime();
		addTransTime(totalTime-temptime);
		totalTime-=beginTime;

		printf("\nafter FFT r(t)=p(t)q(t)\n");
		printf("r(t)=");
		print(r,wlength-1);
		printf("\nuse processor size=%d\n",size);
		printf("total running time=%f(s)\n",totalTime);
		printf("distributed data time=%f(s)\n",transTime);
		printf("parallel compute time=%f(s)\n",totalTime-transTime);
	}

	MPI_Finalize();
}