如何使用SSE或GLSL优化"u [0]*v [0] + u [2]*v [2]"代码行

Question

我有以下功能(来自开源项目"重新导航"):

/// Derives the dot product of two vectors on the xz-plane. (@p u . @p v)
///  @param[in]     u       A vector [(x, y, z)]
///  @param[in]     v       A vector [(x, y, z)]
/// @return The dot product on the xz-plane.
///
/// The vectors are projected onto the xz-plane, so the y-values are ignored.
inline float dtVdot2D(const float* u, const float* v)
{
    return u[0]*v[0] + u[2]*v[2];
}

我通过VS2010 CPU性能测试运行它,它告诉我,在所有重铸代码库中,此功能的代码行u[0]*v[0] + u[2]*v[2]是最热的CPU.

如何进行CPU优化(例如通过SSE或GLSL,如GLM(如果它更容易或更快,在这种情况下更合适)这一行？

编辑:调用显示的上下文:

bool dtClosestHeightPointTriangle(const float* p, const float* a, const float* b, const float* c, float& h) {
    float v0[3], v1[3], v2[3];
    dtVsub(v0, c,a);
    dtVsub(v1, b,a);
    dtVsub(v2, p,a);

    const float dot00 = dtVdot2D(v0, v0);
    const float dot01 = dtVdot2D(v0, v1);
    const float dot02 = dtVdot2D(v0, v2);
    const float dot11 = dtVdot2D(v1, v1);
    const float dot12 = dtVdot2D(v1, v2);

    // Compute barycentric coordinates
    const float invDenom = 1.0f / (dot00 * dot11 - dot01 * dot01);
    const float u = (dot11 * dot02 - dot01 * dot12) * invDenom;
    const float v = (dot00 * dot12 - dot01 * dot02) * invDenom;

Answer 1

在纸上尝试了几件事后,我想出了一些可能对你有用的东西.它是SSE中函数的正确并行/矢量化实现.

然而,它确实需要数据重组,因为我们将同时对4个三角形进行并行处理.

我将逐步分解它并在这里和那里使用指令名称,但请使用C内在函数(_mm_load_ps(),_ mm_sub_ps()等,它们在VC中的xmmintrin.h中) - 当我说到寄存器时只是意味着__m128.

步骤1.

我们根本不需要Y坐标,因此我们设置指向X和Z对的指针.每次调用至少提供4对(即总共4个三角形).我将每个X和Z对称为顶点.

第2步.

使用MOVAPS(要求指针对齐到16位)将每个指针指向的前两个顶点加载到寄存器中.

从a加载的寄存器将如下所示:[a0.x,a0.z,a1.x,a1.z]

第3步.

现在,使用单个减法指令,您可以一次计算2个顶点的增量(v0,v1,v2).

计算v0,v1和v2不仅适用于前2个三角形,还适用于后2个三角形!正如我所说,每个输入应该提供总共4个顶点或8个浮点数.只需重复该数据的步骤2和3.

现在我们有2对vx寄存器,每对包含2个三角形的结果.我将它们称为vx_0(第一对)和vx_1(第二对),其中x为0到2.

第4步.

点产品.为了并行化重心计算(稍后),我们需要4个三角形中的每一个的每个点积的结果,在1个单个寄存器中.

所以在你计算dot01的地方,我们会做同样的事,但一次只能有4个三角形.每个v -register包含2个向量的结果,因此我们首先将它们相乘.

让我们说你的点积函数中的u和v - 参数现在是v0_0和v1_0(计算dot01):

将u和v相乘得到:[(v0_0.x0)*(v1_0.x0),(v0_0.z0)*(v1_0.z0),(v0_0.x1)*(v1_0.x1),(v0_0.z1)*(v1_0.z1)]

由于.x0/.z0和.x1/.z1,这可能看起来很混乱,但是看看在步骤2中加载了什么:a0,a1.

如果现在这仍然感觉模糊,拿起一张纸并从头开始写.

接下来,仍然针对相同的点积,对v0_1和v1_1(即第二对三角形)进行乘法运算.

现在我们有2个寄存器,每个寄存器有2个X和Z对(总共4个顶点),相乘并准备加在一起形成4个独立的点积.SSE3有一个指令来完成这个,它被称为HADDPS:

xmm0 = [A,B,C,D] xmm1 = [E,F,G,H]

HADDPS xmm0,xmm1这样做:

xmm0 = [A + B,C + D,E + F,G + H]

它将从我们的第一个寄存器中取出X和Z对,从第二个寄存器中取出它们,将它们加在一起并将它们存储在目标寄存器的第一,第二,第三和第四个组件中.Ergo:此时我们已经为所有4个三角形提供了这个特殊的点积!

**现在对所有点积重复此过程:dot00等.**

第5步.

最后的计算(据我提供的代码所见)是重心的东西.这是代码中的100%标量计算.但是你的输入现在不是标量点积结果(即单浮点数),它们是SSE向量/寄存器,每个4个三角形都有一个点积.

因此,如果你通过使用对所有4个浮点运算的并行SSE运算来展平矢量化,那么最终会得到1个寄存器(或结果),其中包含4个高度,每个三角形1个.

因为我的午休时间已经过去了,所以我不打算拼出这个,但考虑到我给出的设置/想法,这是最后一步,不应该很难弄明白.

我知道这个想法有点延伸,需要一些爱好来自它上面的代码,也许需要一些铅笔和纸张的质量时间,但它会很快(如果你愿意,你甚至可以添加OpenMP) ).

祝好运 :)

(原谅我的模糊解释,如果需要,我可以掀起功能=))

UPDATE

我已经编写了一个实现,并没有像我预期的那样,主要是因为Y组件涉及到你最初粘贴在你的问题中的代码片段(我查了一下).我在这里所做的不仅是要求你重新排列XZ对中的所有点并按4进给它们,而且还为4个三角形中的每个三角形提供3个指针(对于点A,B和C).从当地的角度来看,这是最快的.我仍然可以修改它以要求来自被叫方的更少侵入性更改,但请让我知道什么是可取的.

然后是一个免责声明:这个代码很简单就像地狱一样(我发现在SSE方面与编译器一起工作得很好......他们可以重新组织,因为看起来合适,x86/x64 CPU也可以分享它们).也是命名,这不是我的风格,它不是任何人的,只是用你认为合适的方式做.

希望它有所帮助,如果不是,我会很乐意再次讨论它.如果这是一个商业项目,那么我也可以选择让我加入;)

无论如何,我把它放在pastebin上:http://pastebin.com/20u8fMEb