這是我們Adreno?工程師Vladislav Shimanskiy 撰寫的Adreno GPU 矩陣乘法系列文章的第二部分，也是最后一個(gè)部分。上一個(gè)部分Vladislav Shimanskiy解釋了Adreno 4xx和5xx GPU系列設(shè)備端矩陣乘法（MM）內(nèi)核函數(shù)和主機(jī)端參考代碼的優(yōu)化實(shí)現(xiàn)相關(guān)概念。本文中，他將結(jié)合代碼分析，詳細(xì)介紹基于OpenCL的主機(jī)代碼和內(nèi)核函數(shù)的實(shí)現(xiàn)。

Vlad Shimanskiy是Qualcomm? GPU計(jì)算解決方案團(tuán)隊(duì)的高級(jí)工程師。

正如我上次在討論問(wèn)題“GPU矩陣乘法存在哪些困難？”時(shí)提到的，由于近來(lái)依賴于卷積的深度學(xué)習(xí)引起廣泛關(guān)注，矩陣乘法（MM）運(yùn)算也在GPU上變得流行起來(lái)。像Adreno GPU這樣的并行計(jì)算處理器是加速此類運(yùn)算的理想選擇。然而，MM算法需要在各個(gè)計(jì)算工作項(xiàng)之間共享大量數(shù)據(jù)。因此，優(yōu)化Adreno的MM算法需要我們利用GPU內(nèi)存子系統(tǒng)。

在OpenCL中實(shí)現(xiàn)

前面已經(jīng)給大家介紹了常用的四種優(yōu)化技術(shù)，這里，我們進(jìn)一步介紹在OpenCL中實(shí)現(xiàn)這些優(yōu)化技術(shù)的主機(jī)參考代碼和內(nèi)核函數(shù)，這些參考代碼和內(nèi)核函數(shù)你將可以直接應(yīng)用到你自己的代碼中。

主機(jī)代碼

首先，我們運(yùn)行防止內(nèi)存復(fù)制的主機(jī)代碼。如前文所述，一個(gè)矩陣通過(guò)TP/L1加載，另一個(gè)矩陣通過(guò)常規(guī)全局內(nèi)存訪問(wèn)路徑加載。

兩個(gè)輸入矩陣中的一個(gè)矩陣用圖像表示方法進(jìn)行表示，即示例代碼中的矩陣B，通過(guò)圖像對(duì)矩陣進(jìn)行抽象，并利用圖像讀取原函數(shù)訪問(wèn)，如第一部分中的圖3所示。對(duì)于其他矩陣，都使用全局內(nèi)存緩沖區(qū)進(jìn)行存儲(chǔ)和訪問(wèn)。這也是為什么為矩陣A和矩陣B應(yīng)用不同的內(nèi)存分配方式的原因。而在矩陣C的訪問(wèn)和表示中，因?yàn)橹恍枰仃?/span>C是寫入數(shù)據(jù)，并且每個(gè)矩陣元素只需要寫一次，到C的流量非常低，所以矩陣C將始終通過(guò)直接路徑訪問(wèn)。

矩陣A和C的內(nèi)存分配

下面例程顯示了如何分配可以通過(guò)直接路徑訪問(wèn)的矩陣A和C，這一點(diǎn)相對(duì)簡(jiǎn)單：

cl::Buffer * buf_ptr = new cl::Buffer(*ctx_ptr, CL_MEM_READ_WRITE | CL_MEM_ALLOC_HOST_PTR, na * ma * sizeof(T));

T * host_ptr = static_cast (queue_ptr->enqueueMapBuffer( *buf_ptr, CL_TRUE, CL_MAP_WRITE, 0, na * ma * sizeof(T)));

lda = na;

圖4通過(guò)L2緩存加載的矩陣的內(nèi)存分配（A和C）

根據(jù)前面介紹，為矩陣A和C分配內(nèi)存中，我們是想得到一個(gè)可以被CPU運(yùn)算訪問(wèn)的主機(jī)指針（CPU指針），并且希望可以通過(guò)該指針對(duì)CPU上的緩沖區(qū)進(jìn)行寫入和讀取操作。因此，上述代碼的第1行中調(diào)用OpenCL的Buffer函數(shù)實(shí)現(xiàn)了內(nèi)存分配，并得到了指向CL緩沖區(qū)的指針。

·???????? 該驅(qū)動(dòng)程序分配一個(gè)緩沖區(qū)。

·???????? CL_MEM_ALLOC_HOST_PTR宏表示該內(nèi)存可以被主機(jī)訪問(wèn)。

·???????? 通過(guò)na和ma我們可以指定矩陣的水平和垂直維度。

注意，這里的內(nèi)存不能使用malloc()函數(shù)在主機(jī)CPU上分配；必須在GPU空間中進(jìn)行分配，并在CPU代碼可以寫入之前，將分配得到的內(nèi)存顯式映射到具有CL API映射函數(shù)的CPU地址空間。

在調(diào)用buffer函數(shù)完成了緩沖區(qū)內(nèi)存分配之后，我們必須得到host_ptr指針，在CPU上通過(guò)該指針可以訪問(wèn)分配的矩陣內(nèi)存。

為了得到host_ptr指針，在圖4所示代碼的第2行中，我們調(diào)用了OpenCL API中的enqueueMapBuffer，使用第1行代碼中得到的緩沖區(qū)指針buf_ptr來(lái)獲得host_ptr指針。enqueueMapBuffer函數(shù)返的host_ptr指針是一個(gè)T類型的指針（示例中T是浮點(diǎn)數(shù)），使用host_ptr指針可以在CPU上對(duì)分配得到的矩陣緩存區(qū)內(nèi)存進(jìn)行讀寫。如果我們已經(jīng)分配了矩陣A，這就是我們用來(lái)傳遞該矩陣的指針。

接著我們看到圖4中代碼的第3行，這里通過(guò)lda 確定矩陣每行使用的內(nèi)存量，以類型T為單位。因此，如果我們?cè)诔绦蛑蟹峙湟粋€(gè)100×100矩陣，則lda將為100個(gè)T類型長(zhǎng)度的內(nèi)存空間。（注意，lda不一定等于矩陣的水平維度；在某些情況下，lda可能與之不同）。

這里，我們?cè)谥鳈C(jī)端將lda、ldb和ldc提交給內(nèi)核，以指定矩陣A、B和C的行距。

矩陣B的內(nèi)存分配（圖像）

接下來(lái)我們來(lái)了解矩陣B是如何分配的，矩陣B的分配比前面介紹的矩陣A和C的分配更復(fù)雜，因?yàn)樵诰仃?/span>B的分配中我們使用了2D圖像。

圖像比緩沖區(qū)限制更加嚴(yán)格。它們通常擁有4個(gè)顏色通道（RGBA），并且在內(nèi)存中為圖像分配內(nèi)存空間的時(shí)候必須保證適當(dāng)?shù)膶?duì)齊。這里，我們先假定一個(gè)圖像，并且圖像的每個(gè)顏色分量是一個(gè)浮點(diǎn)數(shù)。如果我們從矩陣的角度來(lái)觀察圖像，我們希望平展顏色分量。如上所述，為提高效率，我們通過(guò)一個(gè)包括4個(gè)float類型數(shù)據(jù)的向量運(yùn)算來(lái)讀取矩陣，將元素按每4個(gè)float類型打包到圖像像素中。因此，我們?cè)谟?jì)算過(guò)程中必須將矩陣的水平大小除以4，這樣我們表示的才是圖像的像素?cái)?shù)量，具體實(shí)現(xiàn)代碼如下圖5所示：

cl::Image * img_ptr = new cl::Image2D(*ctx_ptr, CL_MEM_READ_WRITE | CL_MEM_ALLOC_HOST_PTR, cl::ImageFormat(CL_RGBA, CL_FLOAT), na/4, ma, 0);

cl::size_t<3> origin;

cl::size_t<3> region;

origin[0] = 0; origin[1] = 0; origin[2] = 0;

region[0] = na/4; region[1] = ma; region[2] = 1;

size_t row_pitch;

size_t slice_pitch;

T * host_ptr = static_cast (queue_ptr->enqueueMapImage( *img_ptr, CL_TRUE, CL_MAP_WRITE, origin, region, &row_pitch, &slice_pitch));

ldb = row_pitch / sizeof(T);

圖5：通過(guò)紋理管道（texture pipe） (B)加載的float32矩陣進(jìn)行內(nèi)存分配

上述代碼中，第1行通過(guò)調(diào)用OpenCL中的Image2D函數(shù)來(lái)分配內(nèi)存，與A和C的內(nèi)存分配一樣，使用了CL_MEM_ALLOC_HOST_PTR宏來(lái)指定分配的內(nèi)存可以從主機(jī)端訪問(wèn)。

分配得到圖像可以從主機(jī)端訪問(wèn)的圖像內(nèi)存后，接著看第8行，通過(guò)enqueueMapImage返回可以在CPU端使用的指針host_ptr（和前面矩陣A和C使用的enqueueMapBuffer類似），并確保我們?cè)?/span>GPU內(nèi)存中分配的圖像區(qū)域?qū)τ?/span>CPU可見。在CPU端可以通過(guò)host_ptr訪問(wèn)到該圖像數(shù)據(jù)。

從CPU調(diào)用內(nèi)核函數(shù)

前面已經(jīng)介紹了如何分配內(nèi)存，接下來(lái)介紹如何從CPU調(diào)用內(nèi)核函數(shù)，該操作包括三個(gè)步驟：

?

·???????? 從CPU中取消映射，使矩陣A和B針對(duì)GPU更新。

·???????? 運(yùn)行內(nèi)核函數(shù)。

·???????? 重新映射，使得矩陣C中的結(jié)果對(duì)于CPU可見。

這個(gè)過(guò)程中我們還必須將A和B的內(nèi)存映射回CPU，以便CPU可以更改這些矩陣；但是，這些更改不能同時(shí)被GPU和CPU獲取，需要一個(gè)同步的過(guò)程。在下面的列表中，我們利用了Snapdragon處理器上的共享虛擬內(nèi)存（SVM）方法來(lái)實(shí)現(xiàn)內(nèi)核函數(shù)運(yùn)行周期和內(nèi)存同步：

// update GPU mapped memory with changes made by CPU

queue_ptr->enqueueUnmapMemObject(*Abuf_ptr, (void *)Ahost_ptr);

queue_ptr->enqueueUnmapMemObject(*Bimg_ptr, (void *)Bhost_ptr);

queue_ptr->enqueueUnmapMemObject(*Cbuf_ptr, (void *)Chost_ptr);

// run kernel

err = queue_ptr->enqueueNDRangeKernel(*sgemm_kernel_ptr, cl::NullRange, global, local, NULL, &mem_event);

mem_event.wait();

// update buffer for CPU reads and following writes

queue_ptr->enqueueMapBuffer( *Cbuf_ptr, CL_TRUE, CL_MAP_READ | CL_MAP_WRITE, 0, m_aligned * n_aligned * sizeof(float));

// prepare mapped buffers for updates on CPU

queue_ptr->enqueueMapBuffer( *Abuf_ptr, CL_TRUE, CL_MAP_WRITE, 0, k_aligned * m_aligned * sizeof(float));

// prepare B image for updates on CPU

cl::size_t<3> origin;

cl::size_t<3> region;

origin[0] = 0; origin[1] = 0; origin[2] = 0;

region[0] = n_aligned/4; region[1] = k_aligned; region[2] = 1;

size_t row_pitch;

size_t slice_pitch;

queue_ptr->enqueueMapImage( *Bimg_ptr, CL_TRUE, CL_MAP_WRITE, origin, region, &row_pitch, &slice_pitch);

圖6：內(nèi)核函數(shù)運(yùn)行周期和內(nèi)存同步過(guò)程

上述代碼實(shí)現(xiàn)分為兩個(gè)部分，其中第一部分是使用enqueueUnmapMemObject函數(shù)調(diào)用取消映射過(guò)程。需要傳遞對(duì)CPU端矩陣做出的所有改變，使其對(duì)于GPU可見，供乘法使用。這是一個(gè)緩存一致性事件：我們分配了矩陣A和B，在CPU端傳播，然后使它們對(duì)GPU可見，而不是復(fù)制內(nèi)存。

完成了第一部分的處理，到了第二部分，GPU現(xiàn)在可以看到分配的矩陣了，并且可以使用。enqueueNDRangeKernel運(yùn)行將對(duì)矩陣進(jìn)行運(yùn)算的內(nèi)核函數(shù)。（經(jīng)驗(yàn)豐富的OpenCL程序員知道如何設(shè)置內(nèi)核函數(shù)的參數(shù)，為簡(jiǎn)潔起見，在此予以省略）。

第二部分的其余部分大同小異，不過(guò)與第一部分相反。內(nèi)核函數(shù)將矩陣乘以矩陣C，因此現(xiàn)在我們需要使矩陣C對(duì)CPU可見。MM運(yùn)算經(jīng)常重復(fù)，因此我們將A和B內(nèi)存映射回CPU，為下一個(gè)運(yùn)算周期做好準(zhǔn)備。在下一次迭代時(shí)，CPU能夠?yàn)?/span>A和B分配新值。

運(yùn)行在GPU上的內(nèi)核函數(shù)代碼

前面已經(jīng)知道了如何進(jìn)行內(nèi)存分配和內(nèi)核函數(shù)的調(diào)用，為了進(jìn)一步了解整個(gè)MM運(yùn)算的性能，我們來(lái)分析運(yùn)行在GPU上的MM運(yùn)算內(nèi)核函數(shù)代碼，這部分代碼說(shuō)明了擁有float 32格式元素的MM運(yùn)算的本質(zhì)。它是BLAS庫(kù)中SGEMM運(yùn)算的簡(jiǎn)化版本，C = αAB + βC，（為簡(jiǎn)潔起見）其中，α= 1和β= 0。

__kernel void sgemm_mult_only(

??? ?????????????????????? __global const float *A,

??? ?????????????????????? const int lda,

??? ?????????????????????? __global float *C,

??? ?????????????????????? const int ldc,

??? ?????????????????????? const int m,

??? ?????????????????????? const int n,

??? ?????????????????????? const int k,

??? ?????????????????????? __read_only image2d_t Bi)

{

??? int gx = get_global_id(0);

??? int gy = get_global_id(1);

if (((gx << 2) < n) && ((gy << 3) < m))

??? {

??????? float4 a[8];

??????? float4 b[4];

??????? float4 c[8];

for (int i = 0; i < 8; i++)

??????? {

??????????? c[i] = 0.0f;

??????? }

int A_y_off = (gy << 3) * lda;

for (int pos = 0; pos < k; pos += 4)

??????? {

??????????? #pragma unroll

??????????? for (int i = 0; i < 4; i++)

??????????? {

??????????????? b[i] = read_imagef(Bi, (int2)(gx, pos + i));

??????????? }

int A_off = A_y_off + pos;

#pragma unroll

??????????? for (int i = 0; i < 8; i++)

??????????? {

??????????????? a[i] = vload4(0, A + A_off);

????????????? ??A_off += lda;

??????????? }

#pragma unroll

??????????? for (int i = 0; i < 8; i++)

??????????? {

??????????????? c[i] += a[i].x * b[0] + a[i].y * b[1] + a[i].z * b[2] + a[i].w * b[3];

??????????? }

}

#pragma unroll

??????? for (int i = 0; i < 8; i++)

??????? {

??????????? int C_offs = ((gy << 3) + i) * ldc + (gx << 2);

??????????? vstore4(c[i], 0, C + C_offs);

??????? }

??? }

}

?

圖7：實(shí)現(xiàn)C = A * B矩陣運(yùn)算的內(nèi)核函數(shù)示例

一般而言，我們會(huì)展開固定大小的循環(huán)，然后將從矩陣A中讀取圖像和數(shù)據(jù)的操作進(jìn)行分組。具體過(guò)程如下：

·???????? 開始時(shí)，我們?cè)O(shè)置了一些限制，確保在處理矩陣時(shí)不致嚴(yán)重限制其維度，因此可以部分占用工作組。每個(gè)工作組水平和垂直地覆蓋一定數(shù)量的micro-tile，但是視乎不同的矩陣維度，我們可能面臨這樣的情況，即macro-tile中的micro-tile僅部分被矩陣占用。因此，我們要跳過(guò)macro-tile未占用部分中的任何運(yùn)算；這就是這個(gè)條件的作用。矩陣維度仍然必須是4x8的倍數(shù)。

·???????? 然后，通過(guò)代碼將矩陣C的元素初始化為零。

·???????? 最外層的for循環(huán)遍歷pos參數(shù)，并包含三個(gè)子循環(huán)：

·???????? 第一個(gè)子循環(huán)中，我們通過(guò)擁有read_imagef函數(shù)的TP/L1讀取矩陣B的元素。

·???????? 第二個(gè)子循環(huán)包含直接從L2讀取的矩陣A的元素值。

·???????? 第三個(gè)子循環(huán)計(jì)算部分點(diǎn)積。

·???????? 注意，為提高效率，所有加載/存儲(chǔ)和ALU操作均使用由4個(gè)float元素構(gòu)成的向量。

通過(guò)上述代碼分析，整個(gè)內(nèi)核函數(shù)可能看起來(lái)比較簡(jiǎn)單，但實(shí)際上它是一個(gè)經(jīng)過(guò)高度優(yōu)化、均衡的運(yùn)算和數(shù)據(jù)大小組合。在使用的過(guò)程中南建議使用-cl-fast-relaxed-math標(biāo)記編譯內(nèi)核函數(shù)。

工作組大小

根據(jù)上述分析，macro-tile是由多個(gè)4×8 micro-tile組成。水平和垂直維度中micro-tile確切數(shù)量由2-D工作組大小確定。通常，最好使用較大的工作組，避免GPU計(jì)算單元利用不足。我們可以使用OpenCL API函數(shù)getWorkGroupInfo查詢最大工作組大小。但是，上邊界為工作組中工作項(xiàng)的總數(shù)。因此，我們?nèi)匀豢梢栽诳偟拇笮〉南拗葡?，自由選擇實(shí)際的維度組成。以下是查找正確大小的一般方法：

·???????? 最小化部分占用工作組的數(shù)量。

·???????? 基于不同大小的矩陣開發(fā)啟發(fā)式算法，并在運(yùn)行時(shí)使用。

·???????? 使用為特殊情況量身定制的內(nèi)核函數(shù)；例如，在矩陣維度特別小的時(shí)候。

·???????? 如果GPU卸載開銷成為瓶頸，就在CPU上完成小型MM運(yùn)算。

開始行動(dòng)

如本文中所示，MM是一項(xiàng)瓶頸運(yùn)算，因此，您需要在OpenCL代碼中利用上述高性能技術(shù)。這是一種加速使用Adreno GPU上內(nèi)存子系統(tǒng)的深度學(xué)習(xí)應(yīng)用的有效方法。

?

更多Qualcomm開發(fā)內(nèi)容請(qǐng)?jiān)斠姡篞ualcomm開發(fā)者社區(qū)。