[譯]在CUDA C/C++中如何衡量程式碼效能

FiBird發表於2019-05-12
C++

本文翻譯自NVIDIA官方部落格Parallel Forall，內容僅供參考，如有疑問請訪問原網站:https://devblogs.nvidia.com/p….

在這個系列的第一篇文章中，我們通過用CUDA C/C++實現SAXPY，學習了CUDA C/C++程式設計的基本要素。在這篇文章中，我們會學習如何衡量這個程式以及其他CUDAC/C++程式的效能。我們在之後的文章中經常用到這種效能度量技術，因為程式的效能優化將會變得越來越重要。

譯者注：這個系列是指原文的系列，並不是筆者的專欄。

CUDA效能度量通常是在主機端進行的，我們既可以使用CPU的計時器也可以使用CUDA專門的計時器。在開始學習效能度量技術之前，我們需要討論一下如何同步主機和裝置之間的操作。

主機-裝置同步

讓我們來看一下上一篇部落格中SAXPY的資料傳輸和核函式啟動的主機端程式碼：

cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

saxpy<<<(N+255)/256, 256>>>(N, 2.0, d_x, d_y);

cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

這裡使用cudaMemcpy進行資料傳輸的方式是同步傳輸(或者是阻塞傳輸)方式。同步資料傳輸直到前面所有釋出的CUDA呼叫全部結束之後才會開始，而且同步資料傳輸結束之後，隨後的CUDA呼叫才會開始。因此上面第三行的saxpy核函式只有到第二行的y到d_y的資料傳輸結束之後才會啟動。而在另一方面，核函式啟動卻是非同步的。一旦核函式被啟動，控制權就立刻返回到CPU，並不會等待核函式執行完成。這樣的話就會對最後一行的裝置到主機資料傳輸產生競態條件(race condition)，但是資料傳輸的阻塞特性會確保核函式執行完成後再開始資料傳輸。

譯者注：這裡的競態條件前面提到過，簡單說就是前面的資料操作還未完成，後面的操作卻又要使用前面的資料，這樣就會導致錯誤的結果。

使用CPU的計時器來計算核函式的執行時間

現在我們來看一下如何使用CPU的計時器來給核函式計時。

cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

t1 = myCPUTimer();
saxpy<<<(N+255)/256, 256>>>(N, 2.0, d_x, d_y);
cudaDeviceSynchronize();
t2 = myCPUTimer();

cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

在上面的程式碼中，我們除了使用一般的主機時間戳函式myCPUTimer()，還用到了顯式的同步障礙 cudaDeviceSynchronize()來阻塞CPU執行，直到裝置上釋出的指令全部執行結束為止。如果沒有這個同步障礙，這個程式碼測試的就是核函式的啟動時間而不是執行時間。

使用CUDA事件計時

使用類似cudaDeviceSynchronize()函式的主機裝置同步點的一個問題就是它會拖延GPU管道(stall GPU pipeline)。基於這個原因，CUDA提供了一個相比CPU計時器更輕量級的選擇，那就是使用CUDA事件API。CUDA事件API包括呼叫事件建立和銷燬函式、事件記錄函式以及以毫秒為單位計算兩個被記錄事件的執行時間的函式。

譯者注：這裡拖延GPU管道(stall GPU pipeline)的直接結果就是造成CPU和GPU輪流執行，而不再是並行執行。於是就使得程式的執行時間等於CPU與GPU時間之和。具體可以參考:https://blogs.msdn.microsoft….

CUDA事件使用的是CUDA streams的概念。一個CUDA流只是一系列在裝置上順序執行的操作。不同流中的操作可以交替執行，在某些情況下甚至可以交疊執行，這個特性可以被用在隱藏主機和裝置間的資料傳輸。(我們會在之後的文章中討論)。到目前為止，我們所有的操作都是在預設的流中進行的，或者0號流(也叫做空流)。

下面的程式碼中，我們使用了CUDA事件API來對SAXPY程式碼進行效能度量。

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

cudaEventRecord(start);
saxpy<<<(N+255)/256, 256>>>(N, 2.0f, d_x, d_y);
cudaEventRecord(stop);

cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

cudaEventSynchronize(stop);
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);

cuda事件是cudaEvent_t型別，通過cudaEventCreate()和cudaEventDestroy()進行事件的建立和銷燬。在上面的程式碼中cudaEventRecord()將事件start和stop放在預設的流中，即0號stream。函式cudaEventSynchronize()用來阻塞CPU執行直到指定的事件被記錄。函式 cudaEventElapsedTime()的第一個引數返回start和stop兩個記錄之間消逝的毫秒時間。這個值的精度大約是0.5ms。

記憶體頻寬

既然我們已經可以精確地測量核函式的執行時間，那麼我們就可以用它來計算頻寬。我們需要使用理論的峰值頻寬和有效記憶體頻寬來評估頻寬效率。

理論頻寬

理論頻寬可以通過產品資料中的硬體規格來計算。例如英偉達Tesla M2050 GPU使用的是時脈頻率為1546MHz視訊記憶體位寬為384-bit的DDR(雙倍資料速率)RAM。

使用這些資料，我們可以計算出英偉達Tesla M2050的理論峰值頻寬是148 GB/sec：

$$BW_{Theoretical}=1546 * 106 * (384/8) * 2 / 109 = 148 GB/s $$

在這個表示式中，我們將記憶體的時脈頻率的單位轉化為Hz，然後乘以視訊記憶體寬度(除以8之後，單位由位元轉化為位元組)，又乘以2是因為該顯示卡的RAM是DDR(雙倍資料速率)。最後我們將結果除以10^9得到以GB/s的計算結果。

有效頻寬

我們是通過計算特定程式的活動時間和程式如何訪問資料來計算機有效頻寬的。我們使用下面的公式：

$$BW_{Effective} = (R_B + W_B) / (t * 109)$$

這裡，$BW_{Effective}$是以GB/s的有效頻寬，$R_B$是每個核函式被讀取的位元組數，$W_B$是每個核函式被寫入的位元組數，$t$是以秒為單位的執行時間。我們可以修改SAXPY例子來計算有效頻寬，下面是完整的程式碼：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>

__global__ void saxpy(int n, float a, float *x, float *y)
{
    int i = blockIdx.x*blockDim.x + threadIdx.x;
    if (i < n) y[i] = a*x[i] + y[i];
}

int main(void)
{
    int N = 20 * (1 << 20);
    float *x, *y, *d_x, *d_y;
    x = (float*)malloc(N*sizeof(float));
    y = (float*)malloc(N*sizeof(float));

    cudaMalloc(&d_x, N*sizeof(float));
    cudaMalloc(&d_y, N*sizeof(float));

    for (int i = 0; i < N; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

    cudaEventRecord(start);

    // Perform SAXPY on 1M elements
    saxpy<<<(N+511)/512, 512>>>(N, 2.0f, d_x, d_y);

    cudaEventRecord(stop);

    cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

    cudaEventSynchronize(stop);
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);

    float maxError = 0.0f;
    for (int i = 0; i < N; i++) {
        maxError = max(maxError, abs(y[i]-4.0f));
    }

    printf("Max error: %f
", maxError);
    printf("Effective Bandwidth (GB/s): %f
", N*4*3/milliseconds/1e6);
}

在上面的頻寬計算(譯者注：即表示式N*4*3/milliseconds/1e6)中，N*4是每次陣列讀或寫的位元組數，因子3的含義是對x的讀以及y的讀和寫共3次讀寫操作。程式執行時間被存在變數milliseconds中，把它作為分母即可算出單位時間的頻寬大小。注意源程式中除了新增了一些計算頻寬的功能外，我們也改變了陣列的大小和塊的大小(譯者注：由於該程式碼來自之前的部落格，所以具體的變化可以對比原來的程式，在這裡)。編譯並執行上面的程式碼，我們可以得到：

$ ./saxpy

Max error: 0.000000

Effective Bandwidth (GB/s): 110.374872

測定計算吞吐量

我們剛剛只演示瞭如何測定頻寬，也叫做資料吞吐量。另一種非常重要的效能指標叫做計算吞度量。一種比較通用的測量計算吞吐量的方法是計算GFLOP/s(Giga-FLoating-point OPerations per second)，代表“每秒10億次的浮點運算數”，這裡的Giga就是千兆，即10^9。對於我們的SAXPY計算，測量有效的吞吐量是很簡單的：每個SAXPY元素都會做一次乘法加法操作，因此是典型的2FLOPS，所以我們可以得到：

$$GFLOP/{s_{Effective}} = 2N / (t * 109)$$

其中，$N$是SAXPY操作的元素個數，$t$是以秒為單位的執行時間。就像理論峰值頻寬一樣，理論峰值$GFLOP/s$也可以從產品資料查到(但是計算它卻很難，因為它具有架構依賴性)。例如，Tesla M2050 GPU的理論單精度浮點峰值吞吐量是$1030GFLOP/s$，而雙精度浮點峰值吞吐量是$515GFLOP/s$。SAXPY每次計算讀取12個位元組，但是僅僅只有一條單獨的乘法加法指令(2 FLOPs)，所以很明顯這(資料吞吐量)就是頻寬限制。而且在這種情況(實際上是大部分情況)下，頻寬是最重要的衡量和優化指標。在更復雜的計算中，FLOPs級別的效能測定是很困難的。因此更普遍的方法是使用分析工具來分析計算吞吐量是否是一個瓶頸。這些應用測出的的常常是問題依賴的吞吐量(而不是架構依賴的)，這其實對使用者會更有用。例如天文學裡每秒百萬次互動作用的N體問題，或者每天納秒級的分子動態模擬。

總結

這篇文章主要介紹瞭如何用CUDA事件API獲取核函式的執行時間。CUDA事件使用GPU計時器，因此避免了與主機裝置同步相關的問題。我們也介紹了有效頻寬和計算吞吐量的效能測定方法，而且也應用這些方法測定了SAXPY例子中核函式的有效頻寬。另外我們也得出，它的記憶體頻寬佔了很大比例，因此在效能測試中，計算有效吞吐量是首要的一步。在之後的文章中，我們會進一步討論在頻寬、指令、或者延遲這些因素中，哪一個是限制程式效能的因素。

CUDA事件也可以用來計算主機和裝置之間資料傳輸的速率，方法很簡單隻要將記錄事件的函式放到cudaMemcpy()呼叫的兩邊就可以了。

如果你在一個很小的GPU上執行文章中的程式碼，那麼如果你沒有減小陣列的大小，你可能會得到一個關於不充足裝置記憶體的錯誤訊息。實際上，我們的例項程式碼目前為止還沒有特別檢查執行時錯誤。在下一篇文章中，我們會學習如何進行錯誤處理以及如何訪問現有裝置來確定已有資源，這樣的話我們就可以寫出更魯棒的程式碼。

相關文章