FFT并行算法与应用-基于MPI(三)

最新推荐文章于 2020-06-11 19:23:15 发布
最新推荐文章于 2020-06-11 19:23:15 发布
阅读量7.9k 收藏 14
点赞数 2
分类专栏： MPI
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Pieces_thinking/article/details/79165292
版权

给出针对一维FFT的并行计算方法，并针对其加以分析推广到二维FFT的计算中。
并行FFT算法：
（1）计算长度为p和n/p2的按位倒置序列；
（2）进行处理机内部的数据交换和外部的数据交换；
（3）在每个处理机中计算长度为s的FFT交换；
（4）将这些长度为s的FFT变换合成为长度为n的FFT变换；

代码：

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

#define MAX_N 50
#define PI    3.1415926535897932
#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->r = c1->r + c2->r;
    result->i = c1->i + c2->i;
}


void comp_multiply(complex_t* result, const complex_t* c1, const complex_t* c2)
{
    result->r = c1->r * c2->r - c1->i * c2->i;
    result->i = c1->r * c2->i + c2->r * c1->i;
}


/*
 * Function:    shuffle
 * Description: 移动f中从beginPos到endPos位置的元素，使之按位置奇偶
 *              重新排列。举例说明:假设数组f，beginPos=2, endPos=5
 *              则shuffle函数的运行结果为f[2..5]重新排列，排列后各个
 *              位置对应的原f的元素为: f[2],f[4],f[3],f[5]
 * Parameters:  f为被操作数组首地址
 *              beginPos, endPos为操作的下标范围
 */
void shuffle(complex_t* f, int beginPos, int endPos)
{
    int i;
    complex_t temp[2*MAX_N];

    for(i = beginPos; i <= endPos; i ++)
    {
        temp[i] = f[i];
    }

    int 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++;
    }
}


/*
 * Function:        evaluate
 * Description: 对复数序列f进行FFT或者IFFT(由x决定)，结果序列为y，
 *          产生leftPos 到 rightPos之间的结果元素
 * Parameters:  f : 原始序列数组首地址
 *          beginPos : 原始序列在数组f中的第一个下标
 *          endPos : 原始序列在数组f中的最后一个下标
 *          x : 存放单位根的数组，其元素为w,w^2,w^3...
 *          y : 输出序列
 *          leftPos : 所负责计算输出的y的片断的起始下标
 *          rightPos : 所负责计算输出的y的片断的终止下标
 *          totalLength : y的长度
 */
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))
    {
        printf("Error in use Polynomial!\n");
        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],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);
        }
    }
}


/*
 * Function:    print
 * Description: 打印数组元素的实部
 * Parameters:  f为待打印数组的首地址
 *              fLength为数组的长度
 */
void print(const complex_t* f, int fLength)
{
    BOOL isPrint = FALSE;
    int i;

    /* f[0] */
    if (abs(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");
}


/*
 * Function:    myprint
 * Description: 完整打印复数数组元素，包括实部和虚部
 * Parameters:  f为待打印数组的首地址
 *              fLength为数组的长度
 */
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");
}


/*
 * Function:   addTransTime
 * Description:累计发送数据所耗费的时间
 * Parameters: toAdd为累加的时间
 */
void addTransTime(double toAdd)
{
    transTime += toAdd;
}


/*
 * Function:    readFromFile
 * Description: 从dataIn.txt读取数据
 */
BOOL readFromFile()
{
    int i;
    FILE* fin = fopen("dataIn.txt", "r");

    if (fin == NULL)
    {
        printf("Cannot find input data file\n"
            "Please create a file \"dataIn.txt\"\n"
            "2\n"
            "1.0  2\n"
            "2.0  -1\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.0;
    }

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

    fclose(fin);

    printf("Read from data file \"dataIn.txt\"\n");
    printf("p(t) = ");
    print(p, variableNum);
    printf("q(t) = ");
    print(q, variableNum);

    return(TRUE);
}


/*
 * Function:    sendOrigData
 * Description: 把原始数据发送给其它进程
 * Parameters:  size为集群中进程的数目
 */
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);
    }
}


/*
 * Function:    recvOrigData
 * Description: 接受原始数据
 */
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;

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

    if(rank == 0)
    {
        /* 0# 进程从文件dataIn.txt读入多项式p,q的阶数和系数序列 */
        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();

        /* 0#进程把多项式的阶数、p、q发送给其它进程 */
        sendOrigData(size);

        /* 累计传输时间 */
        addTransTime(MPI_Wtime() - beginTime);
    }
    else                                          /* 其它进程接收进程0发送来的数据，包括variableNum、数组p和q */
    {
        recvOrigData();
    }
    /* 初始化数组w，用于进行傅立叶变换 */
    int 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);
    }

    /* 划分各个进程的工作范围 startPos ~ stopPos */
    int everageLength = wLength / size;
    int moreLength = wLength % size;
    int startPos = moreLength + rank * everageLength;
    int stopPos  = startPos + everageLength - 1;

    if(rank == 0)
    {
        startPos = 0;
        stopPos  = moreLength+everageLength - 1;
    }

    /* 对p作FFT，输出序列为s，每个进程仅负责计算出序列中 */
    /* 位置为startPos 到 stopPos的元素 */
    evaluate(p, 0, variableNum - 1, w, s, startPos, stopPos, wLength);

    /* 对q作FFT，输出序列为r，每个进程仅负责计算出序列中 */
    /* 位置为startPos 到 stopPos的元素 */
    evaluate(q, 0, variableNum - 1, w, r, startPos, stopPos, wLength);

    /* s和r作点积，结果保存在s中，同样，每个进程只计算自己范围内的部分 */
    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;
    }

    /* 各个进程都把s中自己负责计算出来的部分发送给进程0，并从进程0接收汇总的s */
    if (rank > 0)
    {
        MPI_Send(s + startPos, everageLength * 2, MPI_DOUBLE, 0, S_TAG, MPI_COMM_WORLD);
        MPI_Recv(s, wLength * 2, MPI_DOUBLE, 0, S_TAG2, MPI_COMM_WORLD, &status);
    }
    else
    {
        /* 进程0接收s片断，向其余进程发送完整的s */
        double tempTime = MPI_Wtime();

        for(i = 1; i < size; i ++)
        {
            MPI_Recv(s + moreLength + i * everageLength, everageLength * 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);
    }

    /* swap(w[i],w[(wLength-i)%wLength]) */
    /* 重新设置w，用于作逆傅立叶变换 */
    complex_t temp;
    for(i = 1; i < wLength/2; i ++)
    {
        temp = w[i];
        w[i] = w[wLength - i];
        w[wLength - i] = temp;
    }

    /* 各个进程对s作逆FFT，输出到r的相应部分 */
    evaluate(s, 0, wLength - 1, w, r, startPos, stopPos, wLength);

    /* 各进程把自己负责的部分的r的片断发送到进程0 */
    if (rank > 0)
    {
        MPI_Send(r + startPos, everageLength * 2, MPI_DOUBLE,
            0,R_TAG, MPI_COMM_WORLD);
    }
    else
    {
        /* 进程0接收各个片断得到完整的r，此时r就是两多项式p,q相乘的结果多项式了 */
        double tempTime = MPI_Wtime();

        for(i = 1; i < size; i ++)
        {
            MPI_Recv((r+moreLength+i*everageLength), everageLength * 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 prossor size = %d\n", size);
        printf("Total running time = %f(s)\n", totalTime);
        printf("Distribute data time = %f(s)\n", transTime);
        printf("Parallel compute time = %f(s)\n", totalTime - transTime);
    }
    MPI_Finalize();
}
OpenFPGA
关注
  • 2
    点赞
  • 踩
  • 14
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
专栏目录
FFT算法的并行处理研究
06-10
FFT算法的并行处理研究
fft.zip_fft mpi_fft 并行_mpi并行_并行 fft_并行FFT算法
09-21
基于MPI的并行快速傅里叶变换算法，通过多核并行计算提升计算性能。
全并行结构FFT的FPGA实现
08-03
全并行结构FFT的FPGA实现
利用OpenMP/CUDA/MPI对FFT算法优化程序课程报告
12-17
利用OpenMP/CUDA/MPI三种优化方法，实现对快速傅里叶(FFT)算法优化程序课程报告。。。。。。。。
Fft.rar_fft 并行_fft 并行 mpi_fft编程_并行 fft_并行FFT算法
09-22
fft并行算法实现，C语言开发，用mpi编程完成
基于并行处理的FFT快速算法.pdf
06-06
基于并行处理的FFT快速算法 清晰无遮挡版。
FFT的并行实现
热门推荐
LSEC小陆的博客
12-21 1万+
关键词：快速傅里叶变换 高维FFT 并行计算 快速傅里叶变换简介 离散傅里叶变换 离散傅里叶变换（DFT）一般定义为： Fn≡∑k=0N−1fke−2πink/N F_n \equiv \sum_{k=0}^{N-1}f_ke^{-2\pi ink/N}Fn​≡k=0∑N−1​fk​e−2πink/N 离散傅里叶逆变换可以定义为： fk=1N∑k=0N−1Fne−2πink/N f_...
FFT并行算法与应用-基于MPI(一)
OpenFPGA的博客
01-24 1万+
首先，介绍下FFT算法。 参考： http://www.gatevin.moe/acm/fft%E7%AE%97%E6%B3%95%E5%AD%A6%E4%B9%A0%E7%AC%94%E8%AE%B0/ 博客写的很好，也很基础。 简单介绍以下FFT： 为什么需要FFT 　　第一个问题是为什么要创造FFT，简单的说，为了速度。我们承认DFT很有用，但是我们发现他的速度不是很快，
FFT并行算法与应用-基于MPI(二)
OpenFPGA的博客
01-25 9234
二维串行FFT 　　 二维FFT的是实现方法是先对行做FFT将结果放回该行，然后再对列做FFT结果放在该列，计算完所有的列以后，结果就是响应的二维FFT。 代码： main.c #include "fft12_ifft12.h" #include void main() { int i,j; COMPLEX td[16]={{1,4},{2,3},{3,0},{4,0},{5
CUDA实现FFT并行计算
Liadrinz的博客
06-11 6337
前言   本文是个人学习心得的分享，希望大家在阅读文章后能在评论中一起学习交流！另外还可以访问我的HelloCUDA仓库查看我在学习CUDA中写的一些demo程序。 内容概要 复数的CUDA C++实现 从DFT到FFT FFT蝴蝶操作 CUDA中的分治 FFT的并行化 前置知识 算法基础知识 并行计算基础知识 C++基础知识 CUDA编程基础 离散傅里叶变换的原理 从DFT到FFT   离散傅里叶变化(Discrete Fourier Transform)是傅里叶变换的简化版本，计算机科学家发明出
mpi程序并行算法源码 陈国良书中的全部算法源码
03-30
fft并行算法实现，C语言开发，用mpi编程完成 陈国良并行算法的设计与分析中算法的源码全部
FFT并行MPI实现
04-14
FFT算法的MPI实现源代码，编程风格还行，易于理解
mpi fft mpi fft
07-09
有关mpi的fft，很不错的东西，自己顶一个
并行算法实践－mpi源程序
10-19
含有几十个C语言例子程序 很不错哦 都可以运行
mpi经典代码,包括fft,矩阵相乘等
09-20
mpi经典代码,包括fft,矩阵相乘等，非常适合MPICH学习研究
论文研究-一种在多核嵌入式平台上实现FFT的快速并行算法.pdf
07-22
快速傅里叶变换（fast Fourier transform，FFT）算法是对实时数字信号进行快速分析处理的一种基本方法。针对多核嵌入式实时环境下并行FFT算法进行了研究，以有效提高实时信号处理的速度。提出了一种新的静态多项式FFT算法，充分利用静态多项式奇偶项的不同特点直接代入数据计算，免去了层层迭代的计算过程、减少了运算过程中的通信，提高了并行性能。对算法的理论进行了严密论证，通过嵌入式实时平台上运行测试和仿真实验，证实了在数据分段较短的约束条件下，提出算法较经典的FFT并行算法在时间复杂度上有一定优势。多项式静态FFT算法能够有效提高并行FFT运行速度。
快速傅里叶变换FFT的mpi实现
01-10
一种基于mpi的FFT（快速傅里叶变换）实现 数字信号处理 DSP 离散信号分析
三维湍流Rayleigh-Bénard热对流的高效并行直接求解方法* (2015年)
05-22
采用FFT变换解耦压力泊松方程,将其变换成沿z方向上的块三对角方程组,并利用块三对角方程的MPI与OpenMP联立的大规模高效并行近似解求解方案, 创建了可以高效并行计算的热对流直接求解方法. 通过对该方法并行效率的...
基于matlab实现的空间调制通信过程，包含信号调制、天线选择等发送过程，以及采用最大似然估计的检测过程 .rar
最新发布
05-02
基于matlab实现的空间调制通信过程，包含信号调制、天线选择等发送过程，以及采用最大似然估计的检测过程。.rar
A7-基于FPGA的FFT算法并行优化
05-11
一种常用的FFT并行优化方法是基于蝶形算法（Butterfly Algorithm），即将FFT分解为多个子FFT，每个子FFT通过并行计算来加速整体计算。 具体实现上，可以使用Vivado HLS等工具将FFT算法的C语言代码转换为FPGA可执行...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

OpenFPGA

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值