OpenACC 优化矩阵乘法

▶ 按书上的步骤使用不同的导语优化矩阵乘法

● 已经优化的代码

#include <iostream>
#include <cstdlib>
#include <chrono>

#define SIZE 1024

using namespace std;
using namespace std::chrono;

double a[SIZE][SIZE], b[SIZE][SIZE], c[SIZE][SIZE], d[SIZE][SIZE]; // 四个数组放入 main 里会报错 Segmentation fault (core dumped)

int main()
{
    //int i, j, k;                          // ijk 和 tmp 在循环中使用时才声明会导致运行时间变长
    double tmp;

#pragma acc enter data create(a, b, c)
#pragma acc kernels present(a, b, c)
    {
        for (int i = 0; i < SIZE; i++)      // 初始化 ab
        {
            for (int j = 0; j < SIZE; j++)
                a[i][j] = (double)(i + j);
        }
        for (int i = 0; i < SIZE; i++)      // 初始化 ab
        {
            for (int j = 0; j < SIZE; j++)
                b[i][j] = (double)(i - j);
        }
        for (int i = 0; i < SIZE; i++)      // 每种方法前都要清空 c
        {
            for (int j = 0; j < SIZE; j++)
                c[i][j] = 0.0;
        }
    }

    high_resolution_clock::time_point t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)       // 方法 1，每层循环都 auto
    {        
#pragma acc loop auto                       
        for (int i = 0; i < SIZE; i++)
        {
#pragma acc loop auto
            for (int j = 0; j < SIZE; j++)
            {
#pragma acc loop auto
                for (int k = 0; k < SIZE; k++)
                    c[i][j] += a[i][k] * b[k][j];
            }
        }
    }

    high_resolution_clock::time_point t2 = high_resolution_clock::now();
    duration<double> time = duration_cast<duration<double>>(t2 - t1);    
    printf("Time OpenACC - Auto: %.6lf s.

", time.count());

#pragma acc kernels present(c)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            c[i][j] = 0.0;
    }

    t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)        // 方法 2，外两层 independent，最里层串行
    {
#pragma acc loop independent 
        for (int i = 0; i < SIZE; i++)
        {
#pragma acc loop independent
            for (int j = 0; j < SIZE; j++)
            {
#pragma acc loop independent
                for (int k = 0; k < SIZE; k++)
                    c[i][j] += a[i][k] * b[k][j];
            }
        }
    }

    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenACC - Independent Seq: %.6lf s.

", time.count());

#pragma acc kernels present(c)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            c[i][j] = 0.0;
    }
        
    t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)        // 方法 3，外两层 independent，最里层规约
    {
#pragma acc loop independent
        for (int i = 0; i < SIZE; i++)
        {
#pragma acc loop independent
            for (int j = 0; j < SIZE; j++)
            {
                tmp = 0.0f;
#pragma acc loop reduction(+: tmp)
                for (int k = 0; k < SIZE; k++)
                    tmp += a[i][k] * b[k][j];
                c[i][j] = tmp;
            }
        }
    }

    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenACC - Independent Reduction: %.6lf s.

", time.count());

#pragma acc kernels present(c)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            c[i][j] = 0.0;
    }
    
    t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)        // 方法 4，手动指定 gang 和 vector
    {
#pragma acc loop gang(32)
        for (int i = 0; i < SIZE; i++)
        {
#pragma acc loop vector(16)
            for (int j = 0; j < SIZE; j++)
            {
                tmp = 0.0f;
#pragma acc loop reduction(+: tmp)
                for (int k = 0; k < SIZE; k++)
                    tmp += a[i][k] * b[k][j];
                c[i][j] = tmp;
            }
        }
    }

    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenACC - Gang Vector: %.6lf s.

", time.count());

#pragma acc kernels present(c)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            c[i][j] = 0.0;
    }

    t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)        // 方法 5，分块重排
    {
#pragma acc loop tile(32, 32) 
        for (int i = 0; i < SIZE; i++)
        {
            for (int j = 0; j < SIZE; j++)
            {
                tmp = 0.0f;
#pragma acc loop reduction(+ 
                           : tmp)
                for (int k = 0; k < SIZE; ++k)
                    tmp += a[i][k] * b[k][j];
                c[i][j] = tmp;
            }
        }
    }

    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenACC - tile: %.6lf s.

", time.count());

#pragma acc kernels present(c)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            c[i][j] = 0.0;
    }

    t1 = high_resolution_clock::now();

#pragma acc kernels present(a, b, c)        // 方法 6，合并多层迭代
    {
#pragma acc loop collapse(2) independent
        for (int i = 0; i < SIZE; i++)
        {
            for (int j = 0; j < SIZE; j++)
            {
                tmp = 0.0f;
#pragma acc loop reduction(+: tmp)
                for (int k = 0; k < SIZE; k++)
                    tmp += a[i][k] * b[k][j];
                c[i][j] = tmp;
            }
        }
    }

    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenACC - Collapse: %.6lf s.

", time.count());

#pragma acc exit data copyout(a, b, c)

#pragma omp parallel for shared(d)
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
            d[i][j] = 0.0;
    }

    t1 = high_resolution_clock::now();

#pragma omp parallel for default(none) shared(a, b, d)  // 使用 OpenMP
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            for (int k = 0; k < SIZE; k++)
                d[i][j] += a[i][k] * b[k][j];
        }
    }
    t2 = high_resolution_clock::now();
    time = duration_cast<duration<double>>(t2 - t1);
    printf("Time OpenMP: %.6lf s.

", time.count());

    for (int i = 0; i < SIZE; i++)                      // 检查结果
    {
        for (int j = 0; j < SIZE; j++)
        {
            if (c[i][j] != d[i][j])
                printf("
Error at [%d, %d],c = %f d = %f 
", i, j, c[i][j], d[i][j]);
        }
    }
    return 0;
}

● 输出结果（数据管理优化前）

cuan@CUAN:/mnt/d/Code/OpenACC/OpenACCProject/OpenACCProject$ ./acc.exe

Time OpenACC - Auto: 4.589736 s.

Time OpenACC - Independent Seq: 4.823721 s.

Time OpenACC - Independent Reduction: 3.669336 s.

Time OpenACC - Gang Vector: 3.611391 s.

Time OpenACC - tile: 3.609573 s.

Time OpenACC - Collapse: 3.605792 s.

Time OpenMP: 4.345018 s.

● 输出结果（数据管理优化后）

cuan@CUAN:~/acc$ pgc++ main.cpp -std=c++11 -acc -mp -Minfo -o main.exe
main:
      3, include "chrono"
          31, include "chrono"
              208, Parallel region activated
              212, Parallel region terminated
              217, Parallel region activated
              224, Parallel region terminated
     19, Generating enter data create(b[:][:],c[:][:],a[:][:])
         Generating present(a[:][:],b[:][:],c[:][:])
     20, Loop is parallelizable
     22, Loop is parallelizable
         Generating Tesla code
         20, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
         22, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
     25, Loop is parallelizable
     27, Loop is parallelizable
         Generating Tesla code
         25, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
         27, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
     30, Loop is parallelizable
     32, Loop is parallelizable
         Generating Tesla code
         30, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
         32, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
     32, Memory zero idiom, loop replaced by call to __c_mzero8
     40, Generating present(a[:][:],c[:][:],b[:][:])
     42, Loop is parallelizable
     45, Loop is parallelizable
     48, Complex loop carried dependence of c prevents parallelization
         Loop carried dependence of c prevents parallelization
         Loop carried backward dependence of c prevents vectorization
         Inner sequential loop scheduled on accelerator
         Generating Tesla code
         42, #pragma acc loop gang /* blockIdx.y */
         45, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */
         48, #pragma acc loop seq
     48, Complex loop carried dependence of c prevents parallelization
         Loop carried backward dependence of c prevents vectorization
     56, Generating present(c[:][:])
     59, Loop is parallelizable
     61, Loop is parallelizable
         Generating Tesla code
         59, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
         61, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
     61, Memory zero idiom, loop replaced by call to __c_mzero8
     68, Generating present(a[:][:],c[:][:],b[:][:])
     70, Loop is parallelizable
     73, Loop is parallelizable
     76, Loop is parallelizable
         Generating Tesla code
         70, #pragma acc loop gang /* blockIdx.z */
         73, #pragma acc loop gang, vector(128) /* blockIdx.x threadIdx.x */
         76, #pragma acc loop gang /* blockIdx.y */
     84, Generating present(c[:][:])
     87, Loop is parallelizable
     89, Loop is parallelizable
         Generating Tesla code
         87, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
         89, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
     89, Memory zero idiom, loop replaced by call to __c_mzero8
     96, Generating present(a[:][:],c[:][:],b[:][:])
     98, Loop is parallelizable
    101, Loop is parallelizable
         Generating Tesla code
         98, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
        101, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
        105, #pragma acc loop seq
    101, FMA (fused multiply-add) instruction(s) generated
    105, Loop is parallelizable
    114, Generating present(c[:][:])
    117, Loop is parallelizable
    119, Loop is parallelizable
         Generating Tesla code
        117, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
        119, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
    119, Memory zero idiom, loop replaced by call to __c_mzero8
    126, Generating present(a[:][:],c[:][:],b[:][:])
    128, Loop is parallelizable
    131, Loop is parallelizable
         Generating Tesla code
        128, #pragma acc loop gang(32), vector(8) /* blockIdx.y threadIdx.y */
        131, #pragma acc loop gang, vector(16) /* blockIdx.x threadIdx.x */
        135, #pragma acc loop seq
    135, Loop is parallelizable
    144, Generating present(c[:][:])
    147, Loop is parallelizable
    149, Loop is parallelizable
         Generating Tesla code
        147, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
        149, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
    149, Memory zero idiom, loop replaced by call to __c_mzero8
    156, Generating present(a[:][:],c[:][:],b[:][:])
    158, Loop is parallelizable
    160, Loop is parallelizable
         Generating Tesla code
        158, #pragma acc loop gang, vector tile(32,32) /* blockIdx.x threadIdx.x */
        160,   /* blockIdx.x threadIdx.x tiled */
        165, #pragma acc loop seq
    165, Loop is parallelizable
    174, Generating present(c[:][:])
    177, Loop is parallelizable
    179, Loop is parallelizable
         Generating Tesla code
        177, #pragma acc loop gang, vector(4) /* blockIdx.y threadIdx.y */
        179, #pragma acc loop gang, vector(32) /* blockIdx.x threadIdx.x */
    179, Memory zero idiom, loop replaced by call to __c_mzero8
    186, Generating present(a[:][:],c[:][:],b[:][:])
    188, Loop is parallelizable
    190, Loop is parallelizable
         Generating Tesla code
        188, #pragma acc loop gang, vector(128) collapse(2) /* blockIdx.x threadIdx.x */
        190,   /* blockIdx.x threadIdx.x collapsed */
        194, #pragma acc loop seq
    194, Loop is parallelizable
    208, Generating exit data copyout(c[:][:],b[:][:],a[:][:])
         Parallel loop activated with static block schedule
    210, Memory zero idiom, loop replaced by call to __c_mzero8
    212, Barrier
    217, Parallel loop activated with static block schedule
         FMA (fused multiply-add) instruction(s) generated
    224, Barrier
cuan@CUAN:~/acc$ ./main.exe 
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=22 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=27 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=32 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=48 device=0 threadid=1 num_gangs=8192 num_workers=1 vector_length=128 grid=8x1024 block=128
Time OpenACC - Auto: 0.018726 s.

launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=61 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=76 device=0 threadid=1 num_gangs=32768 num_workers=1 vector_length=128 grid=8x1024x4 block=128
Time OpenACC - Independent Seq: 0.040719 s.

launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=89 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=101 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
Time OpenACC - Independent Reduction: 0.012491 s.

launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=119 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=131 device=0 threadid=1 num_gangs=2048 num_workers=8 vector_length=16 grid=64x32 block=16x8
Time OpenACC - Gang Vector: 0.012314 s.

launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=149 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=160 device=0 threadid=1 num_gangs=1024 num_workers=1 vector_length=1024 grid=1024 block=1024
Time OpenACC - tile: 0.013609 s.

launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=179 device=0 threadid=1 num_gangs=8192 num_workers=4 vector_length=32 grid=32x256 block=32x4
launch CUDA kernel  file=/home/cuan/acc/main.cpp function=main line=190 device=0 threadid=1 num_gangs=8192 num_workers=1 vector_length=128 grid=8192 block=128
Time OpenACC - Collapse: 0.012676 s.

Time OpenMP: 0.504436 s.


Accelerator Kernel Timing data
/home/cuan/acc/main.cpp
  main  NVIDIA  devicenum=0
    time(us): 112,420
    19: compute region reached 1 time
        22: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=38 max=38 min=38 avg=38
            elapsed time(us): total=317 max=317 min=317 avg=317
        27: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=39 max=39 min=39 avg=39
            elapsed time(us): total=50 max=50 min=50 avg=50
        32: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=39 max=39 min=39 avg=39
            elapsed time(us): total=50 max=50 min=50 avg=50
    19: data region reached 3 times
    40: compute region reached 1 time
        48: kernel launched 1 time
            grid: [8x1024]  block: [128]
             device time(us): total=18,705 max=18,705 min=18,705 avg=18,705
            elapsed time(us): total=18,717 max=18,717 min=18,717 avg=18,717
    40: data region reached 2 times
    56: compute region reached 1 time
        61: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=43 max=43 min=43 avg=43
            elapsed time(us): total=176 max=176 min=176 avg=176
    56: data region reached 2 times
    68: compute region reached 1 time
        76: kernel launched 1 time
            grid: [8x1024x4]  block: [128]
             device time(us): total=40,585 max=40,585 min=40,585 avg=40,585
            elapsed time(us): total=40,709 max=40,709 min=40,709 avg=40,709
    68: data region reached 2 times
    84: compute region reached 1 time
        89: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=39 max=39 min=39 avg=39
            elapsed time(us): total=71 max=71 min=71 avg=71
    84: data region reached 2 times
    96: compute region reached 1 time
        101: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=12,456 max=12,456 min=12,456 avg=12,456
            elapsed time(us): total=12,467 max=12,467 min=12,467 avg=12,467
    96: data region reached 2 times
    114: compute region reached 1 time
        119: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=37 max=37 min=37 avg=37
            elapsed time(us): total=63 max=63 min=63 avg=63
    114: data region reached 2 times
    126: compute region reached 1 time
        131: kernel launched 1 time
            grid: [64x32]  block: [16x8]
             device time(us): total=12,295 max=12,295 min=12,295 avg=12,295
            elapsed time(us): total=12,306 max=12,306 min=12,306 avg=12,306
    126: data region reached 2 times
    144: compute region reached 1 time
        149: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=42 max=42 min=42 avg=42
            elapsed time(us): total=187 max=187 min=187 avg=187
    144: data region reached 2 times
    156: compute region reached 1 time
        160: kernel launched 1 time
            grid: [1024]  block: [1024]
             device time(us): total=13,447 max=13,447 min=13,447 avg=13,447
            elapsed time(us): total=13,599 max=13,599 min=13,599 avg=13,599
    156: data region reached 2 times
    174: compute region reached 1 time
        179: kernel launched 1 time
            grid: [32x256]  block: [32x4]
             device time(us): total=41 max=41 min=41 avg=41
            elapsed time(us): total=173 max=173 min=173 avg=173
    174: data region reached 2 times
    186: compute region reached 1 time
        190: kernel launched 1 time
            grid: [8192]  block: [128]
             device time(us): total=12,651 max=12,651 min=12,651 avg=12,651
            elapsed time(us): total=12,669 max=12,669 min=12,669 avg=12,669
    186: data region reached 2 times
    208: data region reached 1 time
        208: data copyout transfers: 3
             device time(us): total=1,963 max=674 min=644 avg=654