2

我正在创建一个简单的矩阵乘法程序，在英特尔至强融核架构上运行。

在多次尝试自动矢量化之后，为了获得更好的性能，我不得不使用 Intel Intrinsics。

到目前为止，矩阵大小是由源代码中的#define 给出的，但是当我尝试在运行时给出它时，性能会大大降低。

源代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <stddef.h>
#include <chrono>
#include <ctime>
#include <mmintrin.h>
#include <xmmintrin.h>  // SSE
#include <pmmintrin.h>  // SSE2
#include <emmintrin.h>  // SSE3
#include <immintrin.h>
#include <zmmintrin.h>

#define ALIGNMENT 64 
#ifndef SIZE
#define SIZE 960
#endif

#define vZero(c) {(c) = _mm512_setzero_pd();}  

#define start_time() \
    auto start = std::chrono::high_resolution_clock::now();
/** Shows the elapsed time. See start_time for usage*/
#define elapsed_time(STRING) \
    auto elapsed = std::chrono::high_resolution_clock::now() - start; \
    long long microseconds = std::chrono::duration_cast<std::chrono::microseconds>(elapsed).count(); \
    printf(#STRING":%lld\n", microseconds);

void recTranspose(double *__restrict__ a, double *__restrict__ aT, const int n, const int k, const int lda, const int ldat){
    if (n*k <= 128) {
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < k; j++) {
                aT[j*ldat+i] = a[i*lda+j];
            }
        }
        //printf("Reached _|_");
        return;
    }
    if(k > n) {
        recTranspose(a, aT, n, (k+1)/2, lda, ldat);
        recTranspose(&a[(k+1)/2], &aT[(k+1)/2*ldat], n, k-((k+1)/2), lda, ldat);
    } else {
        recTranspose(a, aT, (n+1)/2, k, lda, ldat);
        recTranspose(&a[(n+1)/2*lda], &aT[(n+1)/2], n- (n+1)/2, k, lda, ldat);
    }

}
/** Calculates 8 cols and 30 rows of c.*/
inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c, const int size) {
    __m512d c0, c1, c2, c3, c4, c5, c6, c7, c8, c9;
    __m512d c10, c11, c12, c13, c14, c15, c16, c17, c18, c19;
    __m512d c20, c21, c22, c23, c24, c25, c26, c27, c28, c29;


    vZero(c0);  vZero(c1);  vZero(c2);  vZero(c3);  vZero(c4);  vZero(c5);
    vZero(c6);  vZero(c7);  vZero(c8);  vZero(c9);  vZero(c10); vZero(c11);
    vZero(c12); vZero(c13); vZero(c14); vZero(c15); vZero(c16); vZero(c17);
    vZero(c18); vZero(c19); vZero(c20); vZero(c21); vZero(c22); vZero(c23);
    vZero(c24); vZero(c25); vZero(c26); vZero(c27); vZero(c28); vZero(c29);

    __assume_aligned(a, ALIGNMENT);
    __assume_aligned(b, ALIGNMENT);
    __assume_aligned(c, ALIGNMENT);
    __assume(size%16==0);
    for(int i = 0; i < size; i++) {
        const __m512d bv = _mm512_load_pd(b+i*size);
        c0 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+0, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c0);
        c1 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+1, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c1);
        c2 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+2, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c2);
        c3 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+3, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c3);
        c4 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+4, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c4);
        c5 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+5, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c5);
        c6 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+6, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c6);
        c7 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+7, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c7);
        c8 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+8, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c8);
        c9 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+9, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c9);
        c10 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+10, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0),bv, c10);
        c11 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+11, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0),bv, c11);
        c12 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+12, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c12);
        c13 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+13, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c13);
        c14 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+14, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c14);
        c15 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+15, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c15);
        c16 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+16, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c16);
        c17 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+17, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c17);
        c18 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+18, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c18);
        c19 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+19, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c19);
        c20 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+20, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c20);
        c21 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+21, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c21);
        c22 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+22, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c22);
        c23 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+23, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c23);
        c24 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+24, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c24);
        c25 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+25, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c25);
        c26 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+26, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c26);
        c27 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+27, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c27);
        c28 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+28, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c28);
        c29 = _mm512_fmadd_pd(_mm512_extload_pd(a+i*size+29, _MM_UPCONV_PD_NONE, _MM_BROADCAST_1X8, 0), bv, c29);
    }

    _mm512_storenr_pd(c+0*size, c0);
    _mm512_storenr_pd(c+1*size, c1);
    _mm512_storenr_pd(c+2*size, c2);
    _mm512_storenr_pd(c+3*size, c3);
    _mm512_storenr_pd(c+4*size, c4);
    _mm512_storenr_pd(c+5*size, c5);
    _mm512_storenr_pd(c+6*size, c6);
    _mm512_storenr_pd(c+7*size, c7);
    _mm512_storenr_pd(c+8*size, c8);
    _mm512_storenr_pd(c+9*size, c9);

    _mm512_storenr_pd(c+10*size, c10);
    _mm512_storenr_pd(c+11*size, c11);
    _mm512_storenr_pd(c+12*size, c12);
    _mm512_storenr_pd(c+13*size, c13);
    _mm512_storenr_pd(c+14*size, c14);
    _mm512_storenr_pd(c+15*size, c15);

    _mm512_storenr_pd(c+16*size, c16);
    _mm512_storenr_pd(c+17*size, c17);
    _mm512_storenr_pd(c+18*size, c18);
    _mm512_storenr_pd(c+19*size, c19);
    _mm512_storenr_pd(c+20*size, c20);
    _mm512_storenr_pd(c+21*size, c21);
    _mm512_storenr_pd(c+22*size, c22);
    _mm512_storenr_pd(c+23*size, c23);

    _mm512_storenr_pd(c+24*size, c24);
    _mm512_storenr_pd(c+25*size, c25);
    _mm512_storenr_pd(c+26*size, c26);
    _mm512_storenr_pd(c+27*size, c27);
    _mm512_storenr_pd(c+28*size, c28);
    _mm512_storenr_pd(c+29*size, c29);
}



int main(int argc, const char ** argv) {
#ifdef SIZES
    const int size = SIZE;

#else
    const int size = atoi(argv[1]);
#endif
    void* p = malloc((sizeof(double)*5*size*size) + ALIGNMENT-1);
    double *__restrict__ a = (double*)(((size_t)p + ALIGNMENT-1) / ALIGNMENT * ALIGNMENT);
    double *__restrict__ aT = (double*) a+size*size;
    double *__restrict__ b = aT+size*size;
    double *__restrict__ c = b+size*size;
    double *__restrict__ d = c+size*size;
    srand(time(NULL));

    for(int i = 0; i < size; i++) {
        for(int j = 0; j < size; j++) {
            a[i*size+j] = (double) (rand()%20);
        }
        for(int j2=0; j2<size; j2++){
            c[i*size+j2] = 0.0;
        }
    }
    for(int i = 0; i < size; i++) {
        for(int j = 0; j < size; j++) {
            b[i*size+j] = (double) (rand()%20);
        }
    }


    start_time();
    recTranspose(a, aT, size, size, size, size);
    for(int i = 0; i < size; i+=30) {
        for(int j = 0; j < size; j+=8) {
            eightbythirty(&aT[i], &b[j], &c[i*size+j], size);
        }
    }
    elapsed_time();
    double gflops = 2.0*size*size*size*1.0e-03/(microseconds);
    printf("Gflops: %f\n", gflops);

    for(int i = 0; i < size; i++) {
        for(int j = 0; j < size; j++) {
            double s = 0;
            for(int u = 0; u < size; u++) {
                s += a[i*size+u] * b[u*size+j];
            }
            d[i*size+j] = s;
        }
    }

    int error = 0;
    for(int i = 0; i < size; i++) {
        for(int j = 0; j < size; j++) {
            if(abs(c[i*size+j] - d[i*size+j]) > 1) {
                printf("Error at %d %d , %f instead of %f\n", i, j, c[i*size+j], d[i*size+j]);
                error++;
                if(error > 16) return 0;
            }
        }
    }

    printf("OK\n");

}

例如，尺寸为 960（目前它仅适用于尺寸为 30*8 的倍数）：

  • 如果我使用给定大小的编译时间进行编译： icc -mmic -O3 -restrict -std=c++11 -DSIZES -DSIZE=960 mmul.cpp -o mmul.o

    经过时间：0.460745s Gflops：3.840458

  • 如果我使用运行时给定大小进行编译： icc -mmic -O3 -restrict -std=c++11 mmul.cpp -o mmul.o

    经过时间：2.204564s Gflops：0.802640

我认为这可能是 icc 无法识别内存访问模式的预取问题。查看生成的 asm 源代码，“编译时”版本中 vprefetch 指令的数量要高得多。

有趣的事实：在编译时版本中检查乘法的正确结果（代码末尾的两个 for 循环，第 178-197 行）要慢得多！

有什么想法吗？我尝试了#pragma loop_count，但它似乎没用，而且手动内部预取似乎也不是很有效。

提前感谢您的任何回答。

问候，卢卡

4

3 回答 3

1

要检查它是否是预取优化结果，您可以尝试编译您定义大小的版本：

icc -mmic -O3 -restrict -std=c++11 -no-opt-prefetch -DSIZES -DSIZE=960 mmul.cpp -o mmul.o

这将禁用软件预取插入。这是一个合理的假设，因为知道数据集大小将允许更有效的预取。

如果您的问题是预取，您应该尝试更激进的编译器选项 -opt-prefetch=4。另一个想法是使用#pragma prefetch。请参阅此处了解更多信息。

您是否也尝试过#pragma loop_count。在这里检查。

您可以使用 vTune 分析您的代码并检查您的 L1 缓存命中率。对于矩阵乘法，它应该非常高。还要确保使用缓存阻塞。

我还建议使用大量迭代来运行您的基准测试。这将涵盖缓存预热。

我这边的另一个建议是使用数学库来完成这项任务。我看到你的代码是单核的，所以这将导致所有内核的性能约为 250 GFLOPS。那将是大约 25-30% 的效率，这不是很好。使用 MKL 会带来更高的效率。

但是，出于教育目的，这个练习是可以的。

于 2014-11-04T17:18:08.890 回答
1

计算机科学的基本定理指出，任何问题都可以通过另一层间接来解决。

这个想法是将您的代码保留为固定大小的循环，并编写代码以调度到正确的固定大小循环。

首先更改eightbythirty为如下所示：

template<int size>
inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c) {

里面有相同的实现。您可以将其放入，namespace details因为它通常不面向用户。

接下来，包装它：

inline void eightbythirty(double *__restrict__ a, double *__restrict__ b, double * __restrict__ c, const int size_divided_by_240) {
  int size=size_divided_by_240;
  switch( size&7 ) {
    case 0: break;
    case 01: eightbythirty<01>(a,b,c); break;
    case 02: eightbythirty<02>(a,b,c); break;
    case 03: eightbythirty<03>(a,b,c); break;
    case 04: eightbythirty<04>(a,b,c); break;
    case 05: eightbythirty<05>(a,b,c); break;
    case 06: eightbythirty<06>(a,b,c); break;
    case 07: eightbythirty<07>(a,b,c); break;
  }
  a+=(size&7)*8*30;
  b+=(size&7)*8*30;
  c+=(size&7)*8*30;
  switch( (size>>3)&7 ) {
    case 0: break;
    case 01: eightbythirty<1*8>(a,b,c); break;
    case 02: eightbythirty<2*8>(a,b,c); break;
    case 03: eightbythirty<3*8>(a,b,c); break;
    case 04: eightbythirty<4*8>(a,b,c); break;
    case 05: eightbythirty<5*8>(a,b,c); break;
    case 06: eightbythirty<6*8>(a,b,c); break;
    case 07: eightbythirty<7*8>(a,b,c); break;
  }
  a += (size&(7<<3))*8*30;
  b += (size&(7<<3))*8*30;
  c += (size&(7<<3))*8*30;
  switch( (size>>6)&7 ) {
    case 0: break;
    case 01: eightbythirty<1*8*8>(a,b,c); break;
    case 02: eightbythirty<2*8*8>(a,b,c); break;
    case 03: eightbythirty<3*8*8>(a,b,c); break;
    case 04: eightbythirty<4*8*8>(a,b,c); break;
    case 05: eightbythirty<5*8*8>(a,b,c); break;
    case 06: eightbythirty<6*8*8>(a,b,c); break;
    case 07: eightbythirty<7*8*8>(a,b,c); break;
    default:
  }
  a += (size&(7<<6))*8*30;
  b += (size&(7<<6))*8*30;
  c += (size&(7<<6))*8*30;
  int steps = size/8/8/8;
  for( int i = 0; i < steps; ++i ) {
    eightbythirty<512>(a+512*i, b+512*i, c+512*i);
  }
}

这会将您的输入大小分成 3 位块。然后它调用固定大小的实现。上面代码中出现了 4 个分支，其中大部分是简单的跳转表，对于小于 512*8*30 的值。对于大于该值的值，事情主要以 512*8*30 的块完成。

7*3+1 = 22 个你的原始函数的实现，每个都有一个常量size，所以编译器可以完全优化它们。

这通常可以通过元编程来完成，但不值得一次性使用。

当我*(8*30)调用<int size>.eightbythirty

于 2014-11-04T19:58:23.023 回答
0

这个问题也发布在https://software.intel.com/en-us/forums/topic/534856上。software.intel.com 上的回复来自不同的专家组。

于 2014-11-05T17:27:50.367 回答