Xeon :填充性能较慢

Question

我实现了一个简单的n×n矩阵乘法，用OpenMp测试c中相同的性能调优。我的初始代码如下：

#pragma omp parallel for shared(a,b,c) private(h,i,j,k)
    for( i = 0; i < n; i++ )
    {
        for( j = 0; j < n; j++)
        {
            for( k = 0; k < n; k++ )
            {
                a[i*n+j]+= b[i*n+k] * c[k*n+j];

编译器切换j和k循环，所以我得到了ikj算法.我想要实现的第一件事是将每一行填充到64字节，以进行被合并的缓存访问。因此，我计算了每一行所需的额外大小。新尺寸为5904 (参考ISBN-9780124104143)，行长度为5900。我的新代码如下：

#pragma omp parallel for shared(a,b,c) private(h,i,j,k)
    for( i = 0; i < n; i++ )
    {
        for( k = 0; k < n; k++ )
        {
            #pragma simd
            #pragma unroll(4)
            for( j = 0; j < n; j++)
            {
                a[i*pn+j]+= b[i*pn+k] * c[k*pn+j];

其中pn是新的填充行长度。我不得不手动更改循环，因为编译器拒绝这样做。在普通的Xeon处理器上运行这段代码，我的性能几乎和以前一样，也许会好一点。这就是我所期望的。但是当我在Xeon上运行代码时，它大约是inital代码的1/10。

经过进一步的编译器研究后，我注意到，内部循环不再展开和矢量化了。因此，我添加了以下#语用：

#pragma simd
#pragma unroll

矢量化工作得很好，但是余数循环没有展开。性能要好得多，但仍然只有普通版本的1/2左右。

下面是普通代码的编译器(-O3)输出：

  LOOP BEGIN at mm_par.c(75,3)
  remark #25444: Loopnest Interchanged: ( 1 2 3 ) --> ( 1 3 2 )

   LOOP BEGIN at mm_par.c(79,5)
   remark #25460: No loop optimizations reported

    LOOP BEGIN at mm_par.c(77,4)
    remark #15301: PERMUTED LOOP WAS VECTORIZED
    LOOP END

    LOOP BEGIN at mm_par.c(77,4)
    <Remainder>
    remark #25436: completely unrolled by 4  
    LOOP END
   LOOP END
  LOOP END

在这里，使用simd和展开语用的填充的输出：

LOOP BEGIN at mm_ali.c(76,3)
remark #25460: No loop optimizations reported

 LOOP BEGIN at mm_ali.c(78,4)
 remark #25460: No loop optimizations reported

  LOOP BEGIN at mm_ali.c(82,10)
     remark #15301: SIMD LOOP WAS VECTORIZED
  LOOP END
 LOOP END
LOOP END

所以展开就被忽略了。有什么办法强迫它吗？我也问自己，这是否是演出不佳的唯一原因。

编辑:无填充的快速矩阵乘法的assambly如下所示：

 vmovapd   c(%r15,%rbx,8), %zmm1                         #81.28 c1
    vprefetche1 2048+a(%r11,%rbx,8)                         #81.6 c5
    vmovapd   64+c(%r15,%rbx,8), %zmm3                      #81.28 c9
    vprefetch0 768+a(%r11,%rbx,8)                           #81.6 c13
    vmovapd   128+c(%r15,%rbx,8), %zmm4                     #81.28 c17
    vprefetch1 2048+c(%r15,%rbx,8)                          #81.28 c21
    vmovapd   192+c(%r15,%rbx,8), %zmm5                     #81.28 c25
    vprefetch0 768+c(%r15,%rbx,8)                           #81.28 c29
    vfmadd213pd a(%r11,%rbx,8), %zmm0, %zmm1                #81.6 c33
    vprefetche1 2112+a(%r11,%rbx,8)                         #81.6 c37
    vfmadd213pd 64+a(%r11,%rbx,8), %zmm0, %zmm3             #81.6 c41
    vprefetch0 832+a(%r11,%rbx,8)                           #81.6 c45
    vfmadd213pd 128+a(%r11,%rbx,8), %zmm0, %zmm4            #81.6 c49
    vprefetch1 2112+c(%r15,%rbx,8)                          #81.28 c53
    vfmadd213pd 192+a(%r11,%rbx,8), %zmm0, %zmm5            #81.6 c57
    vprefetch0 832+c(%r15,%rbx,8)                           #81.28 c61
    vmovaps   %zmm1, a(%r11,%rbx,8)                         #81.6 c65
    vprefetche1 2176+a(%r11,%rbx,8)                         #81.6 c69
    vmovaps   %zmm3, 64+a(%r11,%rbx,8)                      #81.6 c73
    vprefetch0 896+a(%r11,%rbx,8)                           #81.6 c77
    vmovaps   %zmm4, 128+a(%r11,%rbx,8)                     #81.6 c81
    vprefetch1 2176+c(%r15,%rbx,8)                          #81.28 c85
    vmovaps   %zmm5, 192+a(%r11,%rbx,8)                     #81.6 c89
    vprefetch0 896+c(%r15,%rbx,8)                           #81.28 c93
    vprefetche1 2240+a(%r11,%rbx,8)                         #81.6 c97
    vprefetch0 960+a(%r11,%rbx,8)                           #81.6 c101
    vprefetch1 2240+c(%r15,%rbx,8)                          #81.28 c105
    vprefetch0 960+c(%r15,%rbx,8)                           #81.28 c109
    addq      $32, %rbx                                     #77.4 c113
    cmpq      %rsi, %rbx                                    #77.4 c117
    jb        ..B1.51       # Prob 99%                      #77.4 c117

用于填充的慢乘法器如下所示：

vloadunpackld (%rbx), %zmm0                             #83.6 c1
    addl      $32, %r15d                                    #81.10 c1
    vprefetch1 2048+c(%rcx)                                 #83.30 c5
    vloadunpackhd 64(%rbx), %zmm0                           #83.6 c9
    addq      $256, %rbx                                    #81.10 c9
    vprefetch0 512+c(%rcx)                                  #83.30 c13
    vbroadcastsd b(%r12,%r13,8), %zmm2                      #83.18 c17
    vprefetch1 2112+c(%rcx)                                 #83.30 c21
    vfmadd132pd c(%rcx), %zmm0, %zmm2                       #83.6 c25
    vprefetch0 576+c(%rcx)                                  #83.30 c29
    vpackstoreld %zmm2, (%rsi)                              #83.6 c33
    vprefetch1 2176+c(%rcx)                                 #83.30 c37
    vpackstorehd %zmm2, 64(%rsi)                            #83.6 c41
    addq      $256, %rsi                                    #81.10 c41
    vprefetch0 640+c(%rcx)                                  #83.30 c45
    vloadunpackld (%rdi), %zmm3                             #83.6 c49
    vprefetch1 2240+c(%rcx)                                 #83.30 c53
    vloadunpackhd 64(%rdi), %zmm3                           #83.6 c57
    addq      $256, %rdi                                    #81.10 c57
    vprefetch0 704+c(%rcx)                                  #83.30 c61
    vbroadcastsd b(%r12,%r13,8), %zmm4                      #83.18 c65
    vfmadd132pd 64+c(%rcx), %zmm3, %zmm4                    #83.6 c69
    nop                                                     #83.6 c73
    vpackstoreld %zmm4, (%r8)                               #83.6 c77
    vpackstorehd %zmm4, 64(%r8)                             #83.6 c81
    addq      $256, %r8                                     #81.10 c81
    vloadunpackld (%r9), %zmm5                              #83.6 c85
    vloadunpackhd 64(%r9), %zmm5                            #83.6 c89
    addq      $256, %r9                                     #81.10 c89
    vbroadcastsd b(%r12,%r13,8), %zmm6                      #83.18 c93
    vfmadd132pd 128+c(%rcx), %zmm5, %zmm6                   #83.6 c97
    nop                                                     #83.6 c101
    vpackstoreld %zmm6, (%r10)                              #83.6 c105
    vpackstorehd %zmm6, 64(%r10)                            #83.6 c109
    addq      $256, %r10                                    #81.10 c109
    vloadunpackld (%r11), %zmm7                             #83.6 c113
    vloadunpackhd 64(%r11), %zmm7                           #83.6 c117
    addq      $256, %r11                                    #81.10 c117
    vbroadcastsd b(%r12,%r13,8), %zmm8                      #83.18 c121
    vfmadd132pd 192+c(%rcx), %zmm7, %zmm8                   #83.6 c125
    addq      $256, %rcx                                    #81.10 c129
    movb      %al, %al                                      #83.6 c129
    vpackstoreld %zmm8, (%rdx)                              #83.6 c133
    vpackstorehd %zmm8, 64(%rdx)                            #83.6 c137
    addq      $256, %rdx                                    #81.10 c137
    cmpl      $5888, %r15d                                  #81.10 c141
    jb        ..B1.42       # Prob 99%                      #81.10 c141

这是我的解决方案的完整代码。同样，如果我用n交换np，则性能是n的两倍以上。

#include <sys/time.h>
#include <omp.h>

#ifndef max
#define max(a,b) (((a) (b)) ? (a) : (b))
#define min(a,b) (((a) < (b)) ? (a) : (b))
#endif

#define n 5900
#define pn ((((n*sizeof(double))+63)/64)*(64/sizeof(double))) //padding
#define threadNum 144
#define loops 1

double dtime()
{
    double tseconds = 0.0;
    struct timeval mytime;
    gettimeofday(&mytime, (struct timezone*)0);
    tseconds = (double)(mytime.tv_sec + mytime.tv_usec*1.0e-6);
    return tseconds;

}

double a[n*pn] __attribute__((aligned(64)));
double b[n*pn] __attribute__((aligned(64)));
double c[n*pn] __attribute__((aligned(64)));


main(int argc, char **argv){

    int threadNumber, loopNumber;
    if(argc == 3) 
    {
        threadNumber = atoi(argv[1]);
        loopNumber = atoi(argv[2]);
    } else 
    {
        threadNumber = threadNum;
        loopNumber = loops;     
    }


    double tstart, tstop, ttime;
    int i,j,k,h;

    // initialize matrices
    #pragma omp parallel for
    for(i = 0; i < pn*n; i++)
    {
        a[i]=0.0;
        b[i]=2.0;
        c[i]=2.0;
    }

    omp_set_num_threads(threadNumber); 

    tstart = dtime();

    //parallelize via OpenMP on MIC
    for(h = 0; h < loopNumber; h++){
        #pragma omp parallel for shared(a,b,c) private(h,i,j,k)
        for( i = 0; i < n; i++ )
        {
            for( k = 0; k < n; k++ )
            {           
                #pragma omp simd aligned( a, b, c: 64 )
                for( j = 0; j < n; j++)             
                {
                    a[i*pn+j]+= b[i*pn+k] * c[k*pn +j];
                }
            }
        }

    }

    tstop = dtime();
    double elapsed = tstop - tstart;

    double mFlops = ((double)n)*n*n*2.0*loopNumber/elapsed*1.0e-06;

    #pragma omp parallel
    #pragma omp master
    printf("%d %.3f\n", omp_get_num_threads(), mFlops);

}

Answer 1

您给出的代码片段太小，无法进行适当的测试，因此我可以在没有正确测试的情况下建议解决方案，因此它可能会失败。不管怎么说，这里什么都没有。

您已经说过要填充数据以对齐行，但是，AFAICS从未将此信息传输到编译器，因此编译器无法利用该信息。为此，您有各种解决方案，但是由于您已经有了用于并行化的OpenMP指令，因此使用更多的方法进行矢量化似乎是显而易见的选择。

这样你就能得到这样的东西：

#pragma omp parallel for private( i, j, k )
for( i = 0; i < n; i++ ) {
    double *aa = &a[i * pn];
    double *bb = &b[i * pn];
    for( k = 0; k < n; k++ ) {
        double *cc = &c[k * pn];
        #pragma omp simd aligned( aa, bb, cc: 64 )
        for( j = 0; j < n; j++) {
            aa[j] += bb[k] * cc[j];
        }
    }
}

在这里，我假设您的数组是double，所以您可能需要更改--也就是说，它们不是。我不知道它是否对你有用，但这绝对值得一试。