【新特性】- fft1d 算子功能补全 #1004
Comments
|
FFT1D当前重点规模性能如下： |
|
当前预定开发计划：
|
c2r当前预定开发计划：1、实现C2R 1D调优和测试 5.28-6.13 |
|
R2C预定开发计划： |
mluopR2C FFT 算子开发设计方案
本文档为
1 需求分析1.1 算子需求分析该需求分析为框架原生算子实现功能的需求分析，对于框架原生支持但 MLU-OPS™ 当前版本不支持的功能，需要在 example:
1.2 算子功能和应用场景描述R2C FFT: 对一个长度为N的实数数列进行傅里叶变换，输出长度为 N/2+1的复数数列。因为后半部分结果和前半部分是复共轭的关系，所以该算子仅输出前半部分结果。 1.3 算子输入输出参数要求
1.4 算子限制详细描述与框架需求相比，算子尚有哪些功能/模式/范围/数据类型/xxx 不支持。
example:
1.5 验收标准1.5.1 精度验收标准按照精度验收标准的要求明确本算子的精度标准。 例如：本算子属于纯 IO 类算子，验收标准为 diff3=0。 1.5.2 性能验收标准2 算子接口设计2.1 参考接口
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）。
根据FFT算法原理可以将公式可以利用旋转因子的周期性、对称性以及可约性进行多次递归分解，分解成若干较小规模的DFT进行计算。从而将算法复杂度从 如等式（2），是经过一次基-R大小的FFT分解结果。
在R2C离散傅里叶变换中，除了复数离散傅里叶变换所共有的性质外，还具有复共轭定理的性质。其表示输入均为实数序列时，输出结果具有X(k) 与X(N-k) 是互为共轭复数的关系（实部相等，虚部相反）。因此，为了省略后半部分的冗余数据，R2C FFT中，当输入为N个实数的序列时，输出 1）当k = 0时，其输出序列为X(0)、
2）当 3）当N为偶数时，其输出序列中X(N/2)的虚部也为零。若存在k=N/2R，则
以上就是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、资源分配
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 测试用例设计
其他可根据需要进行补充。算子开发完毕后，补充测试报告链接。 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:
4.2 已经过优化的规模说明此项仅填写未在 4.1 中列出的规模，否则填入 4.1. example:
|
mluopC2R FFT 算子开发设计方案
本文档为
1 需求分析1.1 算子需求分析该需求分析为框架原生算子实现功能的需求分析，对于框架原生支持但 MLU-OPS™ 当前版本不支持的功能，需要在 example:
1.2 算子功能和应用场景描述C2R FFT: 对一个长度为N的复数数列进行逆向傅里叶变换，输出长度为 2N-1的实数数列。 1.3 算子输入输出参数要求
1.4 算子限制详细描述与框架需求相比，算子尚有哪些功能/模式/范围/数据类型/xxx 不支持。
example:
1.5 验收标准1.5.1 精度验收标准按照精度验收标准的要求明确本算子的精度标准。 1.5.2 性能验收标准2 算子接口设计2.1 参考接口
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)$$，$$X_3(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})$$~ 以上就是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 测试用例设计
其他可根据需要进行补充。算子开发完毕后，补充测试报告链接。 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:
4.2 已经过优化的规模说明此项仅填写未在 4.1 中列出的规模，否则填入 4.1. example:
|
|
FFT C2C pr: #1045 (comment) |
|
FFT C2R 确定了新的实现方案，6.10-6.13进行1d的调优和测试的计划预计可以按时完成。 |
|
FFT R2C 已初步实现了小基两级分解，并能正确执行计算。 |
|
在R2C FFT实现过程中，由于输出存在共轭对称部分，因此我们只计算前一半部分的数据。在蝶形计算中，每个section我们只计算一般的蝶形，并对蝶形的后半部分输出进行翻转折叠，转为前半部分对称的内容。如下图： 为最后一级同时计算一个section的5个蝶形，一列为一个蝶形的输出结果。由于R2C FFT的对称性质，我们只计算了前面3个蝶形的数量。在存储之前我们需要对每个蝶形的后半输出进行折叠，折叠成下面的形式进行存储。 在使用__bang_mirror的过程中，发现官方文档中有mirror限制条件，height * sizeof(type) must be divisible by 64 on (m)tp_2xx; |
|
目前r2c fft 1d已基本实现，并能进行一些数据的计算。 |
|
通过列主序实现batch在最低维的性能优化 |
|
通过列主序实现batch在最低维的性能优化 分析：cnnl_copy算子在这个规模下的性能很不好，前后cnnl_copy+行主序的实现中前后copy的时间占了95%以上，而列主序绕过了这个调用，所以性能提升很大。比cufft的性能好的原因，可能是cufft也采用先拷贝成连续数据再进行计算的方式。 |
|
上次review的内容已修改完成，可以再次review了 |
|
除 [768, 90, 180] 外，其他case stride都是contiguous，可以添加几个case，stride 修改成内存上不连续的，类似 [768, 90, 180] 这类 dense stride |
已按要求增加stride的其他case。 |
开发计划可参考以下节点：
The text was updated successfully, but these errors were encountered: