PetrelYy commented 4 months ago

开发计划可参考以下节点：

方案撰写，xx.xx~xx.xx
开发自测，xx.xx~xx.xx
提出 PR/MR，xx.xx~xx.xx
review（ 3个赞），xx.xx~xx.xx
maintainer 合入

squidruge commented 4 months ago

FFT1D当前重点规模性能如下：

squidruge commented 3 months ago

当前预定开发计划：

方案撰写，7.10~7.15
开发自测，7.16~8.28
提出 PR/MR，8.28~8.31
review（ 3个赞），9.1~9.10
maintainer 合入

shouhoo commented 3 months ago

c2r当前预定开发计划：

1、实现C2R 1D调优和测试 5.28-6.13 1.1完成C2R first stage 5.28-6.3 1.2 上传C2R 1D实现方案 5.30 1.3 完成C2R other stage 6.3-6.10 1.4 进行调优和测试 6.10-6.13 2、实现C2R 2D调优和测试 6.14-7.8 2.1 实现调用C2R FFT1D算子 6.14-6.22 2.2上传C2R 2D实现方案 6.16 2.3实现列主序C2R FFT 6.22-7.5 2.4 进行调优和测试 7.5-7.8

nike-tinghai commented 3 months ago

R2C预定开发计划： 1、实现R2C 1D调优和测试 5.28-6.13 1.1 完成R2C first stage 5.28-6.3 1.2上传R2C 1D实现方案 5.30 1.3 完成R2C other stage 6.3-6.10 1.4 进行调优和测试 6.10-6.13 2、实现R2C 2D调优和测试 6.14-7.8 2.1 实现调用R2C FFT1D算子 6.14-6.22 2.2上传R2C 2D实现方案 6.16 2.3实现列主序R2C FFT 6.22-7.5 2.4 进行调优和测试 7.5-7.8

nike-tinghai commented 3 months ago

mluopR2C FFT 算子开发设计方案

mluopR2C FFT 文档基本信息

算子名称	mluopR2C FFT
编制人/日期	zhaoxiang/2024-5-28

修改记录

版本号	修订人	修订日期	修订描述
V1.0	zhaoxiang	2024-5-28	首次提交

内容描述

本文档为mluopR2C FFT算子的设计文档，包括需求分析、接口设计、方案设计、性能优化记录。

算子需求 checklist
算子接口描述
功能描述
框架版本 + 对应源码路径
需求对应网络
网络中用到的规模
是否需要支持原位
是否需要支持 stride 机制
框架单元测试阈值指标（可选）

1 需求分析

1.1 算子需求分析

该需求分析为框架原生算子实现功能的需求分析，对于框架原生支持但 MLU-OPS™ 当前版本不支持的功能，需要在1.4算子限制 章节中显式注明。未明确注明不支持的功能，默认 MLU-OPS™ 全部支持。

example:

算子功能简介	R2C FFT
需求来源	PyTorch/TensorFlow/
应用网络	Conformer
输入数据类型	half, float
输入标量参数	标准数据类型
输入 Shape	input1: [batches, n0]; input2: [batches, 4]
输入 Layout	input: ARRAY
输出数据类型	complex_half, complex_float
输出 Shape	complex_half, complex_float
输出 Layout	ARRAY
模式(可选）	无
是否含有 dim/axis 等类似语义的参数且该参数支持负数/其他特殊处理	通过stride语义来支持dim
是否含有 labels/index 等类似语义的参数且该参数支持负数/界外情况/其他特殊处理	无
是否需要支持原位	是
是否需要支持 stride 机制	是
是否需要支持广播	否
0 元素检查是否直接返回	否
其他特殊需求(在线量化，融合，转数提前等，可选)	无
本次开发优先支持的规模/模式	(90, 180, 768)

1.2 算子功能和应用场景描述

R2C FFT: 对一个长度为N的实数数列进行傅里叶变换，输出长度为 N/2+1的复数数列。因为后半部分结果和前半部分是复共轭的关系，所以该算子仅输出前半部分结果。

1.3 算子输入输出参数要求

参数	类型（输入/输出）	支持类型	物理布局	规模限制
handle	输入		/	无
input1_desc	输入		/	无
input1	输入	half, float	ARRAY	无
input2_desc	输入		/	无
input2	输入	half, float	ARRAY	无
mode	输入		/	无
output_desc	输入		/	无
output	输出		ARRAY	无

1.4 算子限制

详细描述与框架需求相比，算子尚有哪些功能/模式/范围/数据类型/xxx 不支持。使用表格或分段列出均可。

注意：凡是没有在此处列出，但最终被框架检测到的算子限制，均会被视为算子 bug。在此处列出的限制，算子内做好防呆。

example:

限制类型	详细说明
数据类型限制	不支持 input 和 output 同时为 half
布局限制	不支持 NCHW 的 Layout, 需要明确感知 layout 的算子才需要说明限制, 对于不需要感知 layout 的算子, 不要加额外的防呆
规模限制	output_size <= 384KB
功能限制	不支持 xxx 模式/不支持 xx 和 xxx 的模式组合
数据范围限制	half 类型下，数据需满足[xx, xx]范围，否则有精度问题
原位限制	不支持原位
stride 限制	不支持 stride 机制
广播限制	xxx 参数不支持广播
xx 限制	xxx

1.5 验收标准

1.5.1 精度验收标准

按照精度验收标准的要求明确本算子的精度标准。

例如：本算子属于纯 IO 类算子，验收标准为 diff3=0。

1.5.2 性能验收标准

见 MLU-OPS™ 性能验收标准。

2 算子接口设计

2.1 参考接口

FFTW

fftw_plan fftw_plan_dft_r2c_1d(int n0, double *in, fftw_complex *out, unsigned flags);

fftw_plan fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, fftw_complex *out, const int *onembed, int ostride, int odist, unsigned flags);

void fftw_execute_dft_r2c(const fftw_plan p, double *in, fftw_complex *out);

void fftw_execute_dft_c2r(const fftw_plan p, fftw_complex *in, double *out);

void fftw_execute_dft(const fftw_plan p, fftw_complex *in, fftw_complex *out);

void fftw_destroy_plan(fftw_plan plan);

2.2 接口设计

1.FFT plan描述符定义，用于描述执行过程所需的信息:
typedef struct mluopFFTStruct *mluopFFTPlan_t;

2.创建FFT plan的接口函数
mluopStatus_t mluopCreateFFTPlan(mluopFFTPlan_t *fft_plan);

3.plan初始化接口，根据输入描述符input_desc、输出描述符output_desc、变换维度rank，变换规模n[], reservespace大小reservespace_size, 以及workspace大小 workspace_size分配所需内存以及提前计算kernel中所需的常数
mluopStatus_t mluopMakeFFTPlanMany(mluopHandle_t handle,
                                 mluopFFTPlan_t fft_plan,
                                 const mluopTensorDescriptor_t input_desc,
                                 const mluopTensorDescriptor_t output_desc,
                                 const int rank,
                                 const int n[],
                                 size_t *reservespace_size,
                                 size_t *workspace_size);
4.给plan绑定reservespace指针
mluopStatus_t mluopSetFFTReserveArea(mluopHandle_t handle,
                                   mluopFFTPlan_t fft_plan,
                                   void *reservespace);
5.执行接口，可利用创建后的plan多次执行相应FFT
mluopStatus_t mluopExecFFT(mluopHandle_t handle,
                         const mluopFFTPlan_t fft_plan,
                         const void *input,
                         const float scale_factor,
                         void *workspace,
                         void *output,
                         int direction);
6.Destroy接口，负责释放plan与reservespace
mluopStatus_t mluopDestroyFFTPlan(mluopFFTPlan_t fft_plan);

3 实现方案设计

3.1 R2C FFT算法介绍

快速傅里叶变换是离散傅里叶变换的快速实现方式，为高性能计算领域最重要的基础算法之一，被广泛应用于信号处理、偏微分方程求解等领域。快速傅里叶变换的使用，使得处理大量数据和复杂算法成为可能，如气象学中，用于分析天气模式和大气数据等。目前复数FFT 算法的实现已经足够成熟，但对于在音频处理、图像处理等需要实数序列进行离散傅里叶变换计算的场景不能完全进行替换。

离散傅里叶变换是傅里叶变换在时域和频域都呈现离散的形式，其定义为公式（1）。

公式（1）

根据FFT算法原理可以将公式可以利用旋转因子的周期性、对称性以及可约性进行多次递归分解，分解成若干较小规模的DFT进行计算。从而将算法复杂度从降低到O(nlogn)。

如等式（2），是经过一次基-R大小的FFT分解结果。

公式（2）

在R2C离散傅里叶变换中，除了复数离散傅里叶变换所共有的性质外，还具有复共轭定理的性质。其表示输入均为实数序列时，输出结果具有X(k) 与X(N-k) 是互为共轭复数的关系（实部相等，虚部相反）。因此，为了省略后半部分的冗余数据，R2C FFT中，当输入为N个实数的序列时，输出个复数，其中输出序列中第零个复数的虚部为0，若N 为偶数，第N/2个复数的虚部也为0。基于以上的这些内容，可将R2C FFT 蝶形计算kernel 分为了以下三种情况：

1）当k = 0时，其输出序列为X(0)、、... 、X( )，其中与、与X( )等均满足互为共轭复数的关系，满足实数傅里叶变换的要求，因此只需要计算N/2及其之前的数据即可。并且其输入序列是，均为虚部为零的实数，同时与输入序列相乘值不变。对此，面对这种输入实数输出为R/2+1个数的情况，我们将进行特殊处理。由于有奇数基和偶数基，因此为了便于统一，当R为奇数时，令R'为R-1，偶数时为R-1。将k=0 带入到上述FFT分解公式即C2C蝶形公式中，可以得到等式（3），该公式被称为为R2C_I蝶形kernel。

公式（3）

2）当时，该部分的输出序列为X(k)、X(k+N/R)、X(k+2N/R)、... 、X(k+(R-1)N/R)。其数据在输出序列中不存在互为共轭复数的情况满足复数FFT中的C2C蝶形情况，因此每个数据都是有效的，但在输出序列中只需要X(N/2+1) 及其之前的数据即可。因此，该部分依旧采用复数FFT的蝶形进行计算，然后将X(N/2) 后面的数据进行虚部取反，转换为前半部分未计算的输出数据。

3）当N为偶数时，其输出序列中X(N/2)的虚部也为零。若存在k=N/2R，则的虚部也为零。但此时等式（2）中旋转因子系数不为特殊值，因此需要对该部分进行单独计算。此时的输出序列中与、等，都是关于X(N/2) 成共轭复数的关系。因此只需要计算前半部分即可，X(N/2) 以后的数据可以省略计算。由于当R为奇数时，输出序列的最后一个数为X(N/2)，R 为偶数时则最后一个数为。因此，当R 为奇数时，令R''=R，偶数时为R-1，将奇偶数据进行统一。同时，由于旋转因子系数

中的k = N/2R，将旋转因子系数放入前面的旋转因子矩阵中进行结合，则等式（2）可以表示成等式（4）。该部分也被称成R2C_II蝶形kernel。

公式(4)

以上就是R2C FFT 算法的基本原理，在面对一个输入为实数序列的离散傅里叶变换，可通过上述三种蝶形kernel 应对实数DFT 递归分解的所有情况，并在省略冗余数据的计算的同时，也不会像其他算法一样产生额外转换操作。

3.2 MLU上面的R2C FFT实现

为了跟适应MLU的硬件环境，本文对此进行了方案的适配适配设计。主要分为了一下两个方面。

1、蝶形块的分解和网络的选择

在host端，先针对输入的规模进行分解. 对于一个大规模的基可做两个步骤的分解，第一步是针对片下蝶形网络调用的大基分解，第二步是为了进一步降低计算复杂度，针对片下蝶形网络调用的小基对大基进行了进一步分解。大基的大小受限于nram的容量，计算时需要全部加载到nram中。本设计方案中将2048-14000范围内的数设定为了大基，小于这个范围的为小基。

由于stockHam蝶形网络适合向量化变成，因此我们采用stockham作为我们的迭代FFT方案。使其能够更好的进行输入和输出。由于最终的输出是N/2+1个结果，因此在网络中，我们每次均只计算每个section的前一半结果，后半部分进行了省略。

2、蝶形的计算

在蝶形计算的过程中，针对MLU的矩阵计算单元，本方案拟采用矩阵乘的方法进行蝶形计算。在执行计算前，会根据分解情况，对不同大小的基，生成dft_matrix,其存储了蝶形旋转因子矩阵的计算数据和矩阵大小。计算时，会提前将旋转因子矩阵从GDRAM中拷贝到NRAM.

蝶形计算时，first stage的输入是用户输入的实数序列，因此蝶形的计算如公式3，一个实数序列乘以一个后面部分省略的旋转因子矩阵。如图1,为了增加计算带宽，实现时会同时对多组数据进行计算。计算过程中，我们将实部和虚部进行了分开计算，因此这里的out_r为in乘以存储实部数据的旋转因子矩阵，out_i为in乘以存储实部数据的旋转因子矩阵。

在算法原理中每一个section都存在一个R2C_I蝶形，若N能被2R整除，还存在一个R2C_II蝶形。在MLU中，为更好的保证数据处理的连续性，不单独计算R2C_I蝶形和R2C_II蝶形，本方案在第一级计算时就将蝶形计算中虚部为0的数均进行了存储，使其成为一个完整的复数，同时也不对R2C_II蝶形的旋转因子系数进行融合。这样，在first stage后的每一级输入数据均为复数，采用C2C蝶形进行处理。与C2C FFT不同的是，计算后的输出矩阵，其后半部分的输出位于当前section的后半部分，因我们只需要前一半的结果，因此，我们将后半部分进行了虚部取反和像前半部分进行对称折叠，如图2。

3、与C2C FFT实现差异

由于输入的为实数，输出为前N/2+1的结果，因此在网络中，每个section的后半部分都进行了省略，直到最后一级section的计算,利用R2C FFT的这一特性能进一步地提升计算效率。同时由于每个section存储地数据大小变了，输入步长变为了当前section间跨度距离。

蝶形计算方面也存在着一定的差异，首先是第一级针对实数单独进行了蝶形kernel的设计，第一级以后的级数，由于后半部分section的数据不用存储，因此对C2C蝶形的后半部分输出进行了折叠和转换为所需数据。

3.4 性能优化设计

1、资源分配

表项	分配策略
NRAM	保存 node 数据，分 ping pong 两部分，以及dftmtrx数据
WRAM	保存matmul的右矩阵in_r和in_i数据
SRAM	举例：暂时存储 IO 的输入输出数据，卷积滤波张量最先 load 在 SM 上然后 broadcast
DRAM(workspace)	存储计算中的临时数据，如大基所需的旋转因子和输入输出等

2、流水设计建议这一部分细致的排布之后，能完全理清楚整个实现细节，写代码仅仅将方案用代码的形式进行描述，无论是测试用例设计，还是功能调试，性能分析优化都会十分方便，以下是几个建议：

对计算部分，建议写用的是哪种指令，比如如果是加的话，使用 pool 还是用 atomic_add 还是 cycle_add，一定要写清楚，这样可以暴露一些由于硬件指令不熟悉而影响性能的情况。

对 IO 部分，建议写清楚 IO pattern，因为不同的 IO pattern 是影响性能的关键因素，比如连续 load 数据，读跳 stride 进行访存，写跳 stride 进行对齐。

对流水部分，建议按照模板，可以看到每个时间片都在做什么事情，硬件提供的各个流的实际操作，目的是可以看到是否有些没有数据依赖的操作可以藏起来。

3.5 可维护性设计

1、bangc 代码中加入必要的 log 信息，比如输入的规模、数据类型、layout 这些，以及如果出错会导致程序 core dump 的变量，比如 IO 指令的 data_size、dim xyz 的值等，这些信息都是有利于快速定位问题；

2、对每一个函数命名变量命名都有充分的注释；

3、避免魔鬼数字，对于确定的数字尽量使用公共宏来替代。

3.6 测试用例设计

算子在网络中用到的规模： xxx
边界 case： xxx

其他可根据需要进行补充。算子开发完毕后，补充测试报告链接。

3.7 算子防呆检查

列出算子需要做的防呆，比如

1、handle, plan, desc指针为空防呆；

2、rank 为 1， 2， 3；

3、输入输出维度防呆；

4、输入输出数据类型防呆，针对r2c, c2r, c2c分别防呆；

5、batch 大小防呆；

6、execution dtype 数据类型防呆。

7、输入输出stride防呆；

4 算子性能/精度问题 & 优化记录

4.1 当前存在问题的规模说明

只需列出在测试过程中发现的性能/精度异常的规模。

example:

提交日期	问题规模	问题描述	是否已修复
2022-8-26	N >= 128	性能问题	是
	H >= 128	性能问题	否
	W >= 128	精度问题	否
2022-8-26	C >= 128	精度问题	否

4.2 已经过优化的规模说明

此项仅填写未在 4.1 中列出的规模，否则填入 4.1.

example:

提交日期	修复规模	修复问题
2022-8-26	C == 1	性能提升
2022-8-26	C 不对齐	精度修复

shouhoo commented 3 months ago

mluopC2R FFT 算子开发设计方案

mluopC2R FFT 文档基本信息

算子名称	mluopc2R FFT
编制人/日期	徐小虎/2024-5-30

修改记录

版本号	修订人	修订日期	修订描述
V1.0	徐小虎	2024-5-30	首次提交

内容描述

本文档为mluopC2R FFT算子的设计文档，包括需求分析、接口设计、方案设计、性能优化记录。

算子需求 checklist
算子接口描述
功能描述
框架版本 + 对应源码路径
需求对应网络
网络中用到的规模
是否需要支持原位
是否需要支持 stride 机制
框架单元测试阈值指标（可选）

1 需求分析

1.1 算子需求分析

该需求分析为框架原生算子实现功能的需求分析，对于框架原生支持但 MLU-OPS™ 当前版本不支持的功能，需要在1.4算子限制 章节中显式注明。未明确注明不支持的功能，默认 MLU-OPS™ 全部支持。

example:

算子功能简介	c2R FFT
需求来源	PyTorch/TensorFlow/
应用网络	Conformer
输入数据类型	complex_half, complex_float
输入标量参数	标准数据类型
输入 Shape	input: [batches, n]
输入 Layout	input: ARRAY
输出数据类型	half, float
输出 Shape	input: [batches, n*2-1]
输出 Layout	ARRAY
模式(可选）	无
是否含有 dim/axis 等类似语义的参数且该参数支持负数/其他特殊处理	通过stride语义来支持dim
是否含有 labels/index 等类似语义的参数且该参数支持负数/界外情况/其他特殊处理	无
是否需要支持原位	是
是否需要支持 stride 机制	是
是否需要支持广播	否
0 元素检查是否直接返回	否
其他特殊需求(在线量化，融合，转数提前等，可选)	无

1.2 算子功能和应用场景描述

C2R FFT: 对一个长度为N的复数数列进行逆向傅里叶变换，输出长度为 2N-1的实数数列。

1.3 算子输入输出参数要求

参数	类型（输入/输出）	支持类型	物理布局	规模限制
handle	输入		/	无
input1_desc	输入		/	无
input1	输入	half, float	ARRAY	无
input2_desc	输入		/	无
input2	输入	half, float	ARRAY	无
mode	输入		/	无
output_desc	输入		/	无
output	输出		ARRAY	无

1.4 算子限制

详细描述与框架需求相比，算子尚有哪些功能/模式/范围/数据类型/xxx 不支持。使用表格或分段列出均可。

注意：凡是没有在此处列出，但最终被框架检测到的算子限制，均会被视为算子 bug。在此处列出的限制，算子内做好防呆。

example:

限制类型	详细说明
数据类型限制	不支持 input 和 output 同时为 half
布局限制	不支持 NCHW 的 Layout, 需要明确感知 layout 的算子才需要说明限制, 对于不需要感知 layout 的算子, 不要加额外的防呆
规模限制	output_size <= 384KB
功能限制	不支持 xxx 模式/不支持 xx 和 xxx 的模式组合
数据范围限制	half 类型下，数据需满足[xx, xx]范围，否则有精度问题
原位限制	不支持原位
stride 限制	不支持 stride 机制
广播限制	不支持广播

1.5 验收标准

1.5.1 精度验收标准

按照精度验收标准的要求明确本算子的精度标准。

1.5.2 性能验收标准

见 MLU-OPS™ 性能验收标准。

2 算子接口设计

2.1 参考接口

FFTW

fftw_plan fftw_plan_dft_r2c_1d(int n0, double *in, fftw_complex *out, unsigned flags);

fftw_plan fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, fftw_complex *out, const int *onembed, int ostride, int odist, unsigned flags);

void fftw_execute_dft_r2c(const fftw_plan p, double *in, fftw_complex *out);

void fftw_execute_dft_c2r(const fftw_plan p, fftw_complex *in, double *out);

void fftw_execute_dft(const fftw_plan p, fftw_complex *in, fftw_complex *out);

void fftw_destroy_plan(fftw_plan plan);

2.2 接口设计

1.FFT plan描述符定义，用于描述执行过程所需的信息:
typedef struct mluopFFTStruct *mluopFFTPlan_t;

2.创建FFT plan的接口函数
mluopStatus_t mluopCreateFFTPlan(mluopFFTPlan_t *fft_plan);

3.plan初始化接口，根据输入描述符input_desc、输出描述符output_desc、变换维度rank，变换规模n[], reservespace大小reservespace_size, 以及workspace大小 workspace_size分配所需内存以及提前计算kernel中所需的常数
mluopStatus_t mluopMakeFFTPlanMany(mluopHandle_t handle,
                                 mluopFFTPlan_t fft_plan,
                                 const mluopTensorDescriptor_t input_desc,
                                 const mluopTensorDescriptor_t output_desc,
                                 const int rank,
                                 const int n[],
                                 size_t *reservespace_size,
                                 size_t *workspace_size);
4.给plan绑定reservespace指针
mluopStatus_t mluopSetFFTReserveArea(mluopHandle_t handle,
                                   mluopFFTPlan_t fft_plan,
                                   void *reservespace);
5.执行接口，可利用创建后的plan多次执行相应FFT
mluopStatus_t mluopExecFFT(mluopHandle_t handle,
                         const mluopFFTPlan_t fft_plan,
                         const void *input,
                         const float scale_factor,
                         void *workspace,
                         void *output,
                         int direction);
6.Destroy接口，负责释放plan与reservespace
mluopStatus_t mluopDestroyFFTPlan(mluopFFTPlan_t fft_plan);

3 实现方案设计

3.1 C2R FFT算法介绍

快速傅里叶变换是离散傅里叶变换的快速实现方式，为高性能计算领域最重要的基础算法之一，被广泛应用于信号处理、偏微分方程求解等领域。快速傅里叶变换的使用，使得处理大量数据和复杂算法成为可能，如气象学中，用于分析天气模式和大气数据等。目前复数FFT 算法的实现已经足够成熟，但对于在音频处理、图像处理等需要实数序列进行离散傅里叶变换计算的场景不能完全进行替换。

离散傅里叶变换是傅里叶变换在时域和频域都呈现离散的形式，其定义为公式$$X(k)=\sum_{n=0}^{N-1}x(n)W_N^{nk}$$。

根据FFT算法原理可以将公式可以利用旋转因子的周期性、对称性以及可约性进行多次递归分解，分解成若干较小规模的DFT进行计算。从而将算法复杂度从$$O(n^2)$$降低到$$O(nlogn)$$。

如公式（2），是经过一次基-R大小的FFT分解结果。

公式（2）

在C2R离散傅里叶变换中，除了复数离散傅里叶变换所共有的性质外，还具有复共轭定理的性质。C2R FFT中，当输入为N个复数的序列时，本质上是有$$2N-1$$个复数作为输入序列，只不过后$$N-1$$个输入序列被省略了，因为它们与第2个到第N个输入序列存在共轭对称的关系，它们实部相等，虚部互为相反数。此外，C2R FFT 蝶形计算kernel 具有以下特点：

第一点是，当$$k=0$$时，$$X_1(k)$$的输入序列中，$$x(r)$$与$$x(N-r)$$，$$x(2r)$$与$$x(N-2r)$$，...，$$x(\frac{N-r}{2r})$$与$$x(\frac{N-r}{2r}+r)$$为共轭复数，同样可采用C2R蝶形，省略后半部分数据的读取。

第二点是，在$$X_2(k)$$，$$X3(k)$$，...，$$X{(r+1)/2}(k)$$离散傅里叶变换的输入序列，与R2C离散傅里叶变换的输出序列类似，后半部分的输入序列可利用对称性的特点取前半部分中与之相对应的数，再将虚部进行取反转为计算所需要的数据。如$$X_2(k)$$中的输入序列为$$x(ir+1)$$，当$$ir+1>\frac{N-1}{2}$$时，取$$x(ir+1)$$的共轭复数$$x(N-ir-1)$$，并对其虚部取反再进行计算。因此，输入序列中，$$x(\frac{N-r}{2r})$$~ $$x(N-1)$$的数据不需要读取和参与计算，能在一定程度上减少访存开销。

以上就是C2R FFT 算法的基本原理，在面对一个输入为复数序列、输出为实数序列的逆傅里叶变换，可利用上述特点，在省略冗余数据的计算的同时，也不会像其他算法一样产生额外转换操作。

3.2 MLU上面的C2R FFT实现

为了跟适应MLU的硬件环境，本文对此进行了方案的适配适配设计。主要分为了一下两个方面。

1、蝶形块的分解和网络的选择

在host端，先针对输入的规模进行分解. 对于一个大规模的基可做两个步骤的分解，第一步是针对片下蝶形网络调用的大基分解，第二步是为了进一步降低计算复杂度，针对片下蝶形网络调用的小基对大基进行了进一步分解。大基的大小受限于nram的容量，计算时需要全部加载到nram中。本设计方案中将2048-14000范围内的数设定为了大基，小于这个范围的为小基。

由于stockHam蝶形网络适合向量化变成，因此我们采用stockham作为我们的迭代FFT方案。使其能够更好的进行输入和输出。

2、蝶形的计算

在蝶形计算的过程中，针对MLU的矩阵计算单元，本方案拟采用矩阵乘的方法进行蝶形计算。在执行计算前，会根据分解情况，对不同大小的基，生成dft_matrix,其存储了蝶形旋转因子矩阵的计算数据和矩阵大小。计算时，会提前将旋转因子矩阵从GDRAM中拷贝到NRAM.

蝶形计算时，first stage的输入是用户输入的复数，因此蝶形的计算是一个复数序列乘以一个后面部分省略的旋转因子矩阵。如图,为了增加计算带宽，实现时会同时对多组数据进行计算。计算过程中，我们将实部和虚部进行了分开计算，因此这里的out_r为in乘以存储实部数据的旋转因子矩阵，out_i为in乘以存储实部数据的旋转因子矩阵。

在算法原理中每一个section都存在单独的C2R蝶形。在MLU中，为更好的保证数据处理的连续性及充分利用矩阵计算单元，不单独计算C2R蝶形。只有在最后一级的计算中，我们会通过C2R蝶形来使输出变为长度为2N-1的实数序列。

3、与C2C FFT实现差异

由于输入的为实数，输出为前N/2+1的结果，因此在网络中，每个section的后半部分都进行了省略，直到最后一级section的计算,利用R2C FFT的这一特性能进一步地提升计算效率。同时由于每个section存储地数据大小变了，输入步长变为了当前section间跨度距离。

蝶形计算方面也存在着一定的差异，首先是第一级针对实数单独进行了蝶形kernel的设计，第一级以后的级数，由于后半部分section的数据不用存储，因此对C2C蝶形的后半部分输出进行了折叠和转换为所需数据。

3.4 性能优化设计

3.5 可维护性设计

3.6 测试用例设计

算子在网络中用到的规模： xxx
边界 case： xxx

其他可根据需要进行补充。算子开发完毕后，补充测试报告链接。

3.7 算子防呆检查

列出算子需要做的防呆，比如

1、handle, plan, desc指针为空防呆；

2、rank 为 1， 2， 3；

3、输入输出维度防呆；

4、输入输出数据类型防呆，针对r2c, c2r, c2c分别防呆；

5、batch 大小防呆；

6、execution dtype 数据类型防呆。

7、输入输出stride防呆；

4 算子性能/精度问题 & 优化记录

4.1 当前存在问题的规模说明

只需列出在测试过程中发现的性能/精度异常的规模。

example:

提交日期	问题规模	问题描述	是否已修复
2022-8-26	N >= 128	性能问题	是
	H >= 128	性能问题	否
	W >= 128	精度问题	否
2022-8-26	C >= 128	精度问题	否

4.2 已经过优化的规模说明

此项仅填写未在 4.1 中列出的规模，否则填入 4.1.

example:

提交日期	修复规模	修复问题
2022-8-26	C == 1	性能提升
2022-8-26	C 不对齐	精度修复

squidruge commented 3 months ago

FFT C2C pr: https://github.com/Cambricon/mlu-ops/pull/1045#issue-2330261459

shouhoo commented 3 months ago

FFT C2R 确定了新的实现方案，6.10-6.13进行1d的调优和测试的计划预计可以按时完成。

nike-tinghai commented 3 months ago

FFT R2C 已初步实现了小基两级分解，并能正确执行计算。

nike-tinghai commented 3 months ago

在R2C FFT实现过程中，由于输出存在共轭对称部分，因此我们只计算前一半部分的数据。在蝶形计算中，每个section我们只计算一般的蝶形，并对蝶形的后半部分输出进行翻转折叠，转为前半部分对称的内容。如下图： Snipaste_2024-06-11_11-17-44

为最后一级同时计算一个section的5个蝶形，一列为一个蝶形的输出结果。由于R2C FFT的对称性质，我们只计算了前面3个蝶形的数量。在存储之前我们需要对每个蝶形的后半输出进行折叠，折叠成下面的形式进行存储。 Snipaste_2024-06-11_11-17-57

在使用__bang_mirror的过程中，发现官方文档中有mirror限制条件，height sizeof(type) must be divisible by 64 on (m)tp_2xx; width sizeof(type) must be divisible by 128 on (m)tp_2xx; 。当float数据处理中，并行蝶形的规模大于16的时候，得16的整数倍才能正确执行翻转折叠。当并行处理多个section，且一个section中的蝶形数量大于16且不为16倍数时，不能正确翻转。请问什么很好的替代方案吗，已考虑过__bang_rotate指令，但面对下面两个section或多个的情况依旧无法处理。 Snipaste_2024-06-11_11-26-52

nike-tinghai commented 2 months ago

目前r2c fft 1d已基本实现，并能进行一些数据的计算。

squidruge commented 2 months ago

通过列主序实现batch在最低维的性能优化

squidruge commented 2 months ago

通过列主序实现batch在最低维的性能优化最后两行分别为makeFFT1dContiguousInput和makeFFT1dContiguousOutput在行主序中的时间占比，以及列主序相对行主序的加速比

分析：cnnl_copy算子在这个规模下的性能很不好，前后cnnl_copy+行主序的实现中前后copy的时间占了95%以上，而列主序绕过了这个调用，所以性能提升很大。比cufft的性能好的原因，可能是cufft也采用先拷贝成连续数据再进行计算的方式。

nike-tinghai commented 2 months ago

上次review的内容已修改完成，可以再次review了

nike-tinghai commented 2 months ago

fft1d_370_manual2.json fft2d_370_manual.json ifft1d_370_manual2.json ifft2d_370_manual.json irfft1d_370_manual2.json rfft1d_370_manual2.json irfft2d_370_manual.json irfft2d_stride_370_manual.json rfft2d_stride_manual.json rfft2d_370_manual.json .json测试文件上传性能测试结果

PetrelYy commented 2 months ago

除 [768, 90, 180] 外，其他case stride都是contiguous，可以添加几个case，stride 修改成内存上不连续的，类似 [768, 90, 180] 这类 dense stride

nike-tinghai commented 1 month ago

除 [768, 90, 180] 外，其他case stride都是contiguous，可以添加几个case，stride 修改成内存上不连续的，类似 [768, 90, 180] 这类 dense stride

已按要求增加stride的其他case。 fft1d_stride_370_manual2.json fft2d_stride_370_manual.json ifft1d_stride_370_manual2.json ifft2d_stride_370_manual.json irfft1d_stride_370_manual2.json irfft2d_stride_370_manual.json rfft1d_stride_370_manual2.json rfft2d_stride_manual.json .json文件上传测试结果上传

Cambricon / mlu-ops

【新特性】- fft1d 算子功能补全 #1004

c2r当前预定开发计划：

mluopR2C FFT 算子开发设计方案

修改记录

内容描述

算子需求 checklist

1 需求分析

1.1 算子需求分析

1.2 算子功能和应用场景描述

1.3 算子输入输出参数要求

1.4 算子限制

1.5 验收标准

1.5.1 精度验收标准

1.5.2 性能验收标准

2 算子接口设计

2.1 参考接口

2.2 接口设计

3 实现方案设计

3.1 R2C FFT算法介绍

3.2 MLU上面的R2C FFT实现

3.4 性能优化设计

3.5 可维护性设计

3.6 测试用例设计

3.7 算子防呆检查

4 算子性能/精度问题 & 优化记录

4.1 当前存在问题的规模说明

4.2 已经过优化的规模说明

mluopC2R FFT 算子开发设计方案

修改记录

内容描述

算子需求 checklist

1 需求分析

1.1 算子需求分析

1.2 算子功能和应用场景描述

1.3 算子输入输出参数要求

1.4 算子限制

1.5 验收标准

1.5.1 精度验收标准

1.5.2 性能验收标准

2 算子接口设计

2.1 参考接口

2.2 接口设计

3 实现方案设计

3.1 C2R FFT算法介绍

3.2 MLU上面的C2R FFT实现

3.4 性能优化设计

3.5 可维护性设计

3.6 测试用例设计

3.7 算子防呆检查

4 算子性能/精度问题 & 优化记录

4.1 当前存在问题的规模说明

4.2 已经过优化的规模说明