MPI实现fft的迭代算法 源于并行计算——结构。算法。编程中伪码
#include "mpi.h"#include <stdio.h>#include <stdlib.h>#include <math.h>#define T 0typedef struct {double real;double img;} com;double PI;//don't use the omega array,not every node needs a whole copy of omega,not efficientstatic inline void cp(com f,com *t);//copy complexstatic inline void add(com a,com b,com* r);static inline void mul(com a,com b,com* r);static inline void sub(com a,com b,com* r);int br(int src,int size);//bit reversevoid show(com a) {printf("%.4f %.4f \n",a.real,a.img);}int main(int argc,char *argv[]) {PI=atan(1)*4;int k,n;//process id and total number of processesMPI_Init(&argc,&argv);MPI_Comm_rank(MPI_COMM_WORLD,&k);MPI_Comm_size(MPI_COMM_WORLD,&n);//mpi communication objMPI_Request isReq;MPI_Request reReq;MPI_Status s;//datacom a;//inputa.real=k+1;a.img=0;com c;//tempcp(a,&c);com w;//omegacom r;//recieveint l=log(n)/log(2);for(int h=l-1;h>=0;--h) {int p=pow(2,h);int q=n/p;if(k%p==k%(2*p)) {MPI_Issend(&c,2,MPI_DOUBLE,k+p,T,MPI_COMM_WORLD,&isReq);MPI_Irecv(&r,2,MPI_DOUBLE,k+p,T,MPI_COMM_WORLD,&reReq);int time=p*(br(k,l)%q);//compute while recieving and sendingw.real=cos(2*PI*time/n);w.img=sin(2*PI*time/n);MPI_Wait(&reReq,&s);mul(r,w,&r);MPI_Wait(&isReq,&s);add(c,r,&c);} else {MPI_Issend(&c,2,MPI_DOUBLE,k-p,T,MPI_COMM_WORLD,&isReq);MPI_Irecv(&r,2,MPI_DOUBLE,k-p,T,MPI_COMM_WORLD,&reReq);int time=p*(br(k-p,l)%q);//compute while recieving and sendingw.real=cos(2*PI*time/n);w.img=sin(2*PI*time/n);MPI_Wait(&reReq,&s);MPI_Wait(&isReq,&s);mul(c,w,&c);//can't modify until sending comes to an endsub(r,c,&c);}MPI_Barrier(MPI_COMM_WORLD);}MPI_Issend(&c,2,MPI_DOUBLE,br(k,l),T,MPI_COMM_WORLD,&isReq);MPI_Recv(&a,2,MPI_DOUBLE,MPI_ANY_SOURCE,T,MPI_COMM_WORLD,&s);printf("b%d:%.4f %.4f \n",k,a.real,a.img);MPI_Wait(&isReq,&s);MPI_Finalize();return 0;}void cp(com f,com *t) {t->real=f.real;t->img=f.img;}void add(com a,com b,com *c) {c->real=a.real+b.real; c->img=a.img+b.img; } void mul(com a,com b,com *c) { c->real=a.real*b.real-a.img*b.img; c->img=a.real*b.img+a.img*b.real; }void sub(com a,com b,com *c) { c->real=a.real-b.real; c->img=a.img-b.img; } int br(int src,int size) {int tmp=src;int des=0;for(int i=size-1;i>=0;--i) {des=((tmp&1)<<i)|des;tmp=tmp>>1;}return des;}没有使用omega数组，而是直接在使用时计算。因为每个节点最多使用log（n）个omega，计算数组需要计算n个。
碟形中上下两个复数的计算和通信顺序有少许差别，因为可能产生读后写错误。