fft_c.txt

fft（快速傅里叶变换）算法的c语言实现

#include "mpi.h" 
#include <complex> 
#include <math.h> 
#include <iostream> 
#pragma comment(lib,"mpi.lib") 

using namespace std; 

const int MAX_N = 50; //最大结点个数 
const int P_TAG = 100; //几个通信时所用的标签 
const int K_TAG = 101; 
const int D_TAG = 102; 
const int D_TAG2 = 102; 
const double PI = 3.141592653; 

int exp(int, int, int); //计算exp(r, i)，比课本上的函数多了一个参数k，代表i的二进制码的位数 
int symr(int, int, int); //求和i仅第r位不同的j，i的第r位为0 
int getr(int, int, int); //取得i的第r位 

int main(int argc, char *argv[]) 
{ 
    complex<double> p[MAX_N]; 
    complex<double> b[MAX_N]; //结果 
    complex<double> q; 
    complex<double> w[MAX_N]; 
    complex<double> d[10]; 
    complex<double> temp; 
    MPI_Status status; //MPI返回状态 
    int rank; //本处理器的编号 
    int size; //总共的处理器的个数 
    int k = 0; //2^k = size 
    double starttime, endtime; 

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

    if(rank == 0) 
    { 
        starttime = MPI_Wtime(); 
    } 

    //计算w(omega)，w[i]里面存放omega的i次方，计算出来以方便后面使用 
    for(int i = 0; i < size; i ++) 
    { 
        w[i]= complex<double> (cos(i*2*PI/size), sin(i*2*PI/size)); 
    } 

    int dtemp = 0; 
    int temp1 = size; 
    k = 0; 
    while(!dtemp) 
    { 
        dtemp = temp1 % 2; 
        temp1 = temp1 / 2; 
        k++; 
    } 
    if(temp1 == 0) 
    { 
        --k; 
    } 
    else 
    { 
        fprintf(stderr, "处理器不合要求，必须是2的n次方个\n"); 
        fflush(stderr); 
        exit(-1); 
    } 

    //0号处理器接受输入，并将第i个数据发送给处理器i 
    if(rank == 0) 
    { 
        for(int i=0; i<size; ++i) 
        { 
            fprintf(stderr, "现在请输入第%d个复数：", i+1); 
            fflush(stderr); 
            double dtemp; 
            cin>>dtemp; 
            p[i].real(dtemp); 
            cin>>dtemp; 
            p[i].imag(dtemp); 
        } 
        d[0] = p[0]; 
        for(int i = 1; i < size; i ++) //将第i个数据发送给处理器i 
        { 
            MPI_Send(&k, 1, MPI_INT, i, K_TAG, MPI_COMM_WORLD); //k 
            MPI_Send(&p[i], 1, MPI_DOUBLE_COMPLEX, i, P_TAG, MPI_COMM_WORLD); //p   
        } 
    } 
    //其余处理器接收0号处理器的广播信息 
    if(rank != 0) 
    { 
        MPI_Recv(&k, 1, MPI_INT, 0, K_TAG, MPI_COMM_WORLD, &status); 
        MPI_Recv(d, 1, MPI_DOUBLE_COMPLEX, 0, P_TAG, MPI_COMM_WORLD, &status); 
    } 

    MPI_Barrier(MPI_COMM_WORLD); //等待所有处理器均接收完自己的数据 

    for(int r=1; r<=k; ++r) 
    { 
        int j = symr(k, r, rank); 
        int t = getr(k, r, rank); 
        if(t == 1) //r位为1，则发送 
        { 
            MPI_Send(&d[r-1], 1, MPI_DOUBLE_COMPLEX, j, D_TAG, MPI_COMM_WORLD); //向与自己对应的处理器发送上次计算结果 
            MPI_Recv(&d[r], 1, MPI_DOUBLE_COMPLEX, j, D_TAG2, MPI_COMM_WORLD, &status); //从与自己对应的处理器接收本次计算结果 
        } 
        else //r位为0，则接收并进行计算 
        { 
            MPI_Recv(&temp, 1, MPI_DOUBLE_COMPLEX, j, D_TAG, MPI_COMM_WORLD, &status); //从与自己对应的处理器接收上次计算结果 
            int x = 0; 
            x = exp(k, r, rank); 
            d[r] = d[r-1] + w[x]*temp; 
            x = exp(k, r, j); 
            temp = d[r-1] + w[x]*temp; 
            MPI_Send(&temp, 1, MPI_DOUBLE_COMPLEX, j, D_TAG2, MPI_COMM_WORLD); //向与自己对应的处理器发送本次计算结果 
        } 
        MPI_Barrier(MPI_COMM_WORLD); //等待所有处理器第r步执行完成 
    } 

    printf("b[%d] = %f + %f i\n", rank, d[k].real(), d[k].imag()); 

    if(rank == 0) 
    { 
        endtime = MPI_Wtime(); 
        printf("计算完成！共耗时 %f 毫秒\n", endtime - starttime); 
    } 
    MPI_Finalize(); 
} 


//计算exp(r, i)，比课本上的函数多了一个参数k，代表i的二进制码的位数 
int exp(int k, int r, int i) 
{ 
    int result = 0; 
    int *binary = new int[k]; //存放i的二进制码的数组 
    int *temp = new int[k]; 
    int j = 0; 

    for(j=0; j<k; ++j) 
    { 
        binary[k - 1 - j] = i % 2; 
        i /= 2; 
    } 

    for(j=0; j<r; ++j) 
    { 
        temp[r - 1 - j] = binary[j]; 
    } 
    for(j=r; j<k; ++j) 
    { 
        temp[j] = 0; 
    } 

    for(j=0; j<k; ++j) 
    { 
        result *= 2; 
        result += temp[j]; 
    } 
    delete [] binary; 
    delete [] temp; 
    return result; 
} 

//求和i仅第r位不同的j，i的第r位为0 
int symr(int k, int r, int i) 
{ 
    int result = 0; 
    int *temp = new int[k]; //存放i的二进制码的数组 
    int j = 0; 

    for(j=0; j<k; ++j) 
    { 
        temp[k - 1 - j] = i % 2; 
        i /= 2; 
    } 

    if(temp[r - 1] == 0) 
    { 
        temp[r - 1] = 1; 
    } 
    else 
    { 
        temp[r - 1] = 0; 
    } 

    for(j=0; j<k; ++j) 
    { 
        result *= 2; 
        result += temp[j]; 
    } 

    delete [] temp; 
    return result; 
} 

//取得i的第r位 
int getr(int k, int r, int i) 
{ 
    int *temp = new int[k]; //存放i的二进制码的数组 

    for(int j=0; j<k; ++j) 
    { 
        temp[k - 1 - j] = i % 2; 
        i /= 2; 
    } 

    int result = temp[r - 1]; 
    delete [] temp; 

    return result; 
}