SSE 内存访问答案

【问题标题】：SSE memory accessSSE 内存访问
【发布时间】：2015-10-28 17:41:43
【问题描述】：

我需要使用 SSE 执行高斯消除，但我不确定如何从 128 位寄存器（每个寄存器存储 4 个元素）访问每个元素（32 位）。这是原始代码（不使用 SSE）：

unsigned int i, j, k;

for (i = 0; i < num_elements; i ++)             /* Copy the contents of the A matrix into the U matrix. */
    for(j = 0; j < num_elements; j++)
        U[num_elements * i + j] = A[num_elements*i + j];


for (k = 0; k < num_elements; k++){             /* Perform Gaussian elimination in place on the U matrix. */
    for (j = (k + 1); j < num_elements; j++){   /* Reduce the current row. */

        if (U[num_elements*k + k] == 0){
            printf("Numerical instability detected. The principal diagonal element is zero. \n");
            return 0;
        }

        /* Division step. */
        U[num_elements * k + j] = (float)(U[num_elements * k + j] / U[num_elements * k + k]);
    }

    U[num_elements * k + k] = 1;             /* Set the principal diagonal entry in U to be 1. */

    for (i = (k+1); i < num_elements; i++){
        for (j = (k+1); j < num_elements; j++)
            /* Elimnation step. */
            U[num_elements * i + j] = U[num_elements * i + j] -\
                                      (U[num_elements * i + k] * U[num_elements * k + j]);

        U[num_elements * i + k] = 0; 
    } 
}

好的，我在这段代码中遇到了分段错误[核心转储]。我是 SSE 的新手。有人可以帮忙吗？谢谢。

 int i,j,k;
 __m128 a_i,b_i,c_i,d_i;

for (i = 0; i < num_rows; i++)
{
for (j = 0; j < num_rows; j += 4)
{
    int index = num_rows * i + j;
   __m128 v = _mm_loadu_ps(&A[index]); // load 4 x floats
   _mm_storeu_ps(&U[index], v);         // store 4 x floats

}
}
for (k = 0; k < num_rows; k++){  

 a_i= _mm_load_ss(&U[num_rows*k+k]);         


    for (j = (4*k + 1); j < num_rows; j+=4){
               b_i= _mm_loadu_ps(&U[num_rows*k+j]);// Reduce the currentrow. 

    if (U[num_rows*k+k] == 0){
    printf("Numerical instability detected.);

        }

        /* Division step. */
        b_i =    _mm_div_ps(b_i, a_i);
  }

    a_i = _mm_set_ss(1);           

    for (i = (k+1); i < num_rows; i++){
  d_i= _mm_load_ss(&U[num_rows*i+k]);
        for (j = (4*k+1); j < num_rows; j+=4){
           c_i= _mm_loadu_ps(&U[num_rows*i+j]); /* Elimnation step. */
        b_i= _mm_loadu_ps(&U[num_rows*k+j]);    
            c_i = _mm_sub_ps(c_i, _mm_mul_ss(b_i,d_i));
        }
       d_i= _mm_set_ss(0); 
    } 
  }

【问题讨论】：

在您最初的尝试中犯了很多基本错误 - 也许先尝试一些更简单的东西，以便学习 SIMD 编码的基础知识？
x86 标签 wiki (stackoverflow.com/tags/x86/info) 有链接。如果您找到一个好的介绍/概述/入门指南，请编辑 wiki 中的链接，或在此处发表评论，我会查看它。大多数链接都是参考资料，当您了解基础知识时很有用，但必须查看哪些指令的确切作用以及可用的内容。
是的，我看到了错误。感谢您的链接。
哦，对不起。是的，矩阵是大小为 num_elements 的正方形。我应该将它用作扁平的二维数组，就像保罗说的那样，目标是学习 SIMD。谢谢。我会尝试以此为基础
@PeterCordes，x86 标签的链接和描述太棒了！这是一个很好的资源。

标签： c sse simd gaussian

【解决方案1】：

为了让你开始，你的第一个循环应该更像这样：

for (i = 0; i < num_elements; i++)
{
    for (j = 0; j < num_elements; j += 4)
    {
        int index = num_elements * i + j;
        __m128i v = _mm_loadu_ps((__m128i *)&A[index]); // load 4 x floats
        _mm_storeu_ps((__m128i *)&U[index], v);         // store 4 x floats
    }
}

这假定num_elements 是4 的倍数，并且A 和U 都没有正确对齐。

【讨论】：

如果他将矩阵存储在单个数组中，为什么要嵌套循环？
@Bob__：为什么不呢？ ;-) 这可能是一个“扁平化”的二维数组，当然他也可以在这种情况下使用memcpy，但由于目标是学习矩阵的 SIMD 编码，因此使用两个嵌套循环更有意义例子。
没关系，但他的原始代码也是错误的。检查循环的边界。不应该是 for (i = 0; i < num_elements; i ++) U[i] = A[i];
@Bob__：我猜矩阵是正方形的，大小为 num_elements x num_elements（除非它们被展平） - 如果 OP 可以澄清（并在必要时修复它们的标量代码）然后我会相应地更新我的答案。