使用 FMA（融合乘法）指令进行复数乘法

Question

我想利用可用的融合乘法加法/减法 CPU 指令来帮助在适当大小的数组上进行复杂的乘法运算。本质上，基本数学如下所示：

void ComplexMultiplyAddToArray(float* pDstR, float* pDstI, const float* pSrc1R, const float* pSrc1I, const float* pSrc2R, const float* pSrc2I, int len)
{
    for (int i = 0; i < len; ++i)
    {
        const float fSrc1R = pSrc1R[i];
        const float fSrc1I = pSrc1I[i];
        const float fSrc2R = pSrc2R[i];
        const float fSrc2I = pSrc2I[i];

        //  Perform complex multiplication on the input and accumulate with the output
        pDstR[i] += fSrc1R*fSrc2R - fSrc1I*fSrc2I;
        pDstI[i] += fSrc1R*fSrc2I + fSrc2R*fSrc1I;
    }
}

您可能会看到，数据是结构化的，其中我们有单独的实数和虚数数组。现在，假设我有以下函数可用作分别执行 ab+c 和 ab-c 的单个指令的内在函数：

float fmadd(float a, float b, float c);
float fmsub(float a, float b, float c);

天真地，我可以看到我可以用一个 fmadd 和一个 fmsub 替换 2 个乘法、一个加法和一个减法，如下所示：

//  Perform complex multiplication on the input and accumulate with the output
pDstR[i] += fmsub(fSrc1R, fSrc2R, fSrc1I*fSrc2I);
pDstI[i] += fmadd(fSrc1R, fSrc2I, fSrc2R*fSrc1I);

这导致了非常适度的性能改进，以及我认为的准确性，但我认为我真的错过了一些可以代数修改数学的东西，这样我就可以替换更多的 mult/add 或 mult/sub组合。在每一行中，都有一个额外的加法和一个额外的乘法，我觉得我可以转换为单个 fma，但令人沮丧的是，如果不更改操作顺序并得到错误的结果，我无法弄清楚如何做到这一点。任何有想法的数学专家？

就这个问题而言，目标平台可能并不那么重要，因为我知道这些指令存在于各种平台上。

Answer 1

我发现以下（有一点帮助）似乎可以得出正确的答案：

pDstR[i] = fmsub(fSrc1R, fSrc2R, fmsub(fSrc1I, fSrc2I, pDstR[i]));
pDstI[i] = fmadd(fSrc1R, fSrc2I, fmadd(fSrc2R, fSrc1I, pDstI[i]));

但奇怪的是，它在 AVX 上的性能提升不如使用半 FMA 保留数学的真实结果部分，而是使用完整 FMA 让假想结果：

pDstR[i] += fmsub(fSrc1R, fSrc2R, fSrc1I*fSrc2I);
pDstI[i] = fmadd(fSrc1R, fSrc2I, fmadd(fSrc2R, fSrc1I, pDstI[i]));

感谢大家的帮助。

Answer 2

这是一个好的开始。您可以再减少一项：

//  Perform complex multiplication on the input and accumulate with the output
pDstR[i] += fmsub(fSrc1R, fSrc2R, fSrc1I*fSrc2I);
pDstI[i] += fmadd(fSrc1R, fSrc2I, fSrc2R*fSrc1I);

这里你可以在虚部的计算中使用另一个fmadd：

pDstI[i] = fmadd(fSrc1R, fSrc2I, fmadd(fSrc2R, fSrc1I, pDstI[i]));

同样你可以对实部做同样的事情，但你需要否定这个论点。如果这使事情变得更快或更慢，很大程度上取决于您正在处理的架构的微时序：

pDstR[i] = fmsub(fSrc1R, fSrc2R, fmadd(fSrc1I, fSrc2I, -pDstR[i]));

顺便说一句，如果您使用restrict 关键字将目标数组声明为非别名，您可能会获得进一步的性能改进。现在编译器必须假设 pDstR 和 pDstI 可能重叠或指向同一块内存。这将阻止编译器在写入 pDstR[i] 之前加载 pDstI[i]。

之后，如果编译器还没有这样做，一些仔细的循环展开也可能会有所帮助。检查编译器的汇编输出！