使用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;
    }
}

您可能会看到，数据是在我们具有独立数量和虚构数字的单独数组的情况下结构的。现在，假设我具有以下功能作为内在的单个指令，分别执行A b c和a b-c：

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

天真地，我可以看到我可以替换2次，一个添加，一个添加和一个fmadd和一个fmsub，例如：

//  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);

我认为这会导致非常适中的性能改进，但我认为我确实缺少可以通过代数修改数学的东西，以便我可以替换几个多/add/add/mult/sub组合。在每一行中，都有一个额外的添加，我觉得我可以将其转换为单个FMA，但是令人沮丧的是，我不知道如何在不更改操作顺序并获得错误的结果的情况下弄清楚如何做。有想法的数学专家吗？

出于问题的目的，目标平台可能并不重要，因为我知道各种平台上都存在这些说明。

Answer 1

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

//  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关键字将目标数组声明为非异差，则您 May 获得进一步的性能改进。现在，编译器必须假设PDSTR和PDSTI可以重叠或指向相同的内存。这将阻止编译器在将写入PDSTR [I]之前加载PDSTI [I]。

之后，如果编译器尚未这样做，则仔细的循环封闭也可能会有所帮助。检查编译器的汇编器输出！

Answer 2

我发现以下内容（有所帮助）似乎会导致正确答案：

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

但奇怪的是，在AVX上的性能并没有像使用半数学的实际结果部分那样提高avx的性能，但是使用完整的FMA，拥有虚构的结果：

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

感谢大家的帮助。