向量化和并行化Xeon Phi

Question

我正在寻找一个简单的例子，在Xeon Phi上使用矢量化和并行化，这比只使用Xeon-Xeon有更好的性能。你能帮我个忙吗？

我正在尝试下一个示例。我注释了run on only-Xeon的第14、18和19行，以及Xeon-Phi的uncoment这些行，但只有-Xeon的性能比Xeon-phi更好

1.void main(){
2.double *a, *b, *c;
3.int i,j,k, ok, n=100;
4.int nPadded = ( n%8 == 0 ? n : n + (8-n%8) );
5.ok = posix_memalign((void**)&a, 64, n*nPadded*sizeof(double));
6.ok = posix_memalign((void**)&b, 64, n*nPadded*sizeof(double));
7.ok = posix_memalign((void**)&c, 64, n*nPadded*sizeof(double));
8.for(i=0; i<n; i++)
9.{
10.    a[i] = (int) rand();
11.    b[i] = (int) rand();
12.    c[i] = 0.0;
13.}
14.#pragma offload target(mic) in(a,b:length(n*nPadded)) inout(c:length(n*nPadded))
15.#pragma omp parallel for
16.for( i = 0; i < n; i++ )
17.    for( k = 0; k < n; k++ )
18.        #pragma vector aligned
19.        #pragma ivdep
20.        for( j = 0; j < n; j++ ){
21.                c[i*nPadded+j] = c[i*nPadded+j] + a[i*nPadded+k]*b[k*nPadded+j]        
22.}

Answer 1

关于自动矢量化的前几个词。自动矢量化的优点是简单。你需要设置一些关键字，然后魔术就会发生，编译器会为你快速编写代码。如果你想走这条路，试试这个manual。

这种方法的缺点是没有简单的方法来理解编译器是如何工作的。在矢量化报告中，你会看到"LOOP In“或"LOOP In NOT VECTORIZED”。但是，如果你想真正理解你的代码是如何工作的，唯一的方法就是查看你的程序集。这不是获取程序集的问题。你需要用-fcode-asm编译程序。但我认为，如果你需要阅读汇编来检查“简单的自动矢量化”方法是如何工作的，那就没有那么简单了。

自动矢量化的替代方案是内部函数(实际上，这不是单一的替代方案)。想想内部函数，比如用C函数包装的汇编。许多内部函数内部包装了单个程序集命令。

我推荐使用这个intrinsics guide。

所以我的简单方法是：

1. 实现单线程引用。您将使用它检查intrinsics version.
2. Implement SSE intrinsics版本的正确性。上交所的内部功能要简单得多，可以在至强Phi的Xeon.
3. Implement AVX-512版本上进行测试。
4. 测量您的速度。

让我们用你的程序来做吧。你的程序有很多不同之处：

1. I使用浮点数代替posix_memalign.
2. I。
3. I使用_mm_malloc代替double。假设n除以16没有余数(在AVX-512矢量寄存器中是16个浮点数)。在这个例子中我没有使用循环剥离。
4. 我使用本机模式而不是卸载模式。KNL是可引导的，所以没有必要使用卸载模式anymore.
5. Also，我认为你的程序是不正确的，因为它在一个时刻修改了来自多个线程的c数组。但让我们认为这并不重要，我们只需要一些计算工作。

我的代码工作时间：

英特尔至强5680

• 参考计算时间: 97.677505秒
• 内部计算时间: 6.189296秒

英特尔至强Phi (KNC) SE10X

• 参考计算时间: 199.0 seconds
• Intrinsics计算时间:2.78秒

代码：

#include <stdio.h>
#include <omp.h>
#include <math.h>
#include "immintrin.h"
#include <assert.h>

#define F_E_Q(X,Y,N) (round((X) * pow(10, N)-(Y) * pow(10, N)) == 0)

void reference(float* a, float* b, float* c, int n, int nPadded);
void intrinsics(float* a, float* b, float* c, int n, int nPadded);

char *test(){
  int n=4800;
  int nPadded = n;

  assert(n%16 == 0);

  float* a = (float*) _mm_malloc(sizeof(float)*n*nPadded, 64);
  float* b = (float*) _mm_malloc(sizeof(float)*n*nPadded, 64);
  float* cRef = (float*) _mm_malloc(sizeof(float)*n*nPadded, 64);
  float* c = (float*) _mm_malloc(sizeof(float)*n*nPadded, 64);
  assert(a != NULL);
  assert(b != NULL);
  assert(cRef != NULL);
  assert(c != NULL);

  for(int i=0, max = n*nPadded; i<max; i++){
    a[i] = (int) rand() / 1804289408.0;
    b[i] = (int) rand() / 1804289408.0;
    cRef[i] = 0.0;
    c[i] = 0.0;
  }
  debug_arr("a", "%f", a, 0, 9, 1);
  debug_arr("b", "%f", b, 0, 9, 1);
  debug_arr("cRef", "%f", cRef, 0, 9, 1);
  debug_arr("c", "%f", c, 0, 9, 1);

  double t1 = omp_get_wtime();
  reference(a, b, cRef, n, nPadded);
  double t2 = omp_get_wtime();
  debug("reference calc time: %f", t2-t1);

  t1 = omp_get_wtime();
  intrinsics(a, b, c, n, nPadded);
  t2 = omp_get_wtime();
  debug("Intrinsics calc time: %f", t2-t1);

  debug_arr("cRef", "%f", cRef, 0, 9, 1);
  debug_arr("c", "%f", c, 0, 9, 1);

  for(int i=0, max = n*nPadded; i<max; i++){
    assert(F_E_Q(cRef[i], c[i], 2));
  }

  _mm_free(a);                
  _mm_free(b);
  _mm_free(cRef);
  _mm_free(c);
  return NULL;
}

void reference(float* a, float* b, float* c, int n, int nPadded){
  for(int i = 0; i < n; i++ )
    for(int k = 0; k < n; k++ )
      for(int j = 0; j < n; j++ )
        c[i*nPadded+j] = c[i*nPadded+j] + a[i*nPadded+k]*b[k*nPadded+j];        
}

#if __MIC__

void intrinsics(float* a, float* b, float* c, int n, int nPadded){
  #pragma omp parallel for
  for(int i = 0; i < n; i++ )
    for(int k = 0; k < n; k++ )
      for(int j = 0; j < n; j+=16 ){
        __m512 aPart = _mm512_extload_ps(a + i*nPadded+k, _MM_UPCONV_PS_NONE, _MM_BROADCAST_1X16, _MM_HINT_NONE);
        __m512 bPart = _mm512_load_ps(b + k*nPadded+j);
        __m512 cPart = _mm512_load_ps(c + i*nPadded+j);
        cPart = _mm512_add_ps(cPart, _mm512_mul_ps(aPart, bPart));
        _mm512_store_ps(c + i*nPadded+j, cPart);
      }
}

#else

void intrinsics(float* a, float* b, float* c, int n, int nPadded){
  #pragma omp parallel for
  for(int i = 0; i < n; i++ )
    for(int k = 0; k < n; k++ )
      for(int j = 0; j < n; j+=4 ){
        __m128 aPart = _mm_load_ps1(a + i*nPadded+k);
        __m128 bPart = _mm_load_ps(b + k*nPadded+j);
        __m128 cPart = _mm_load_ps(c + i*nPadded+j);
        cPart = _mm_add_ps(cPart, _mm_mul_ps(aPart, bPart));
        _mm_store_ps(c + i*nPadded+j, cPart);
      }
}

#endif