Edited at

C++ AMPによるGPGPU入門

More than 3 years have passed since last update.

この記事はC++ Advent Calender 2015の3日目の記事です.

前の日はadatcheyさんでした.

環境構築やGPGPU自体の知識に関しては別の記事を立てましたのでそちらを参照してください.

環境構築編

前記事を要約すると,GPGPUとはGPUを数値計算に利用するための技術であり,C++ AMPはCUDAやOpenCLといったGPGPU向けAPIよりもより高レイヤーで,ハードウェアが抽象化された非常に使いやすいライブラリだということです.

さて, C++ AMPに関してですが,覚えるべき最大の事柄はGPUに送るカーネル(GPUに実行してもらうモジュール),データの条件です.


並列処理の基本

以下のコードを見ていただきましょう.

Hello Worldを出力するサンプルコードです.(kalmar sandboxより引用)

#include <iostream>

#include <amp.h>
using namespace concurrency;
int main() {
int v[11] = {'G', 'd', 'k', 'k', 'n', 31, 'v', 'n', 'q', 'k', 'c'};

array_view<int> av(11, v);
parallel_for_each(av.get_extent(), [=](index<1> idx) restrict(amp) {
av[idx] += 1;
});

for(unsigned int i = 0; i < av.get_extent().size(); i++)
std::cout << static_cast<char>(av(i));
return 0;
}

今は深く考えず, parallel_for_each関数にのみ着目してください.

この関数こそ, カーネルをGPUで実行するためのC++ AMPの中心となる関数です.

parallel_for_each関数に渡されるラムダ式は,アクセラレータ(GPUデバイスと考えてよい)で実行されます.

parallel_for_eachの引数は

parallel_for_each(アクセラレータオブジェクト(省略可能),スレッド数,カーネル)

となっています.スレッド数はメモリ上のN次元空間の境界を指定するベクトルであるconcurrecy::extentクラスで指定します.基本は各次元の要素数と同じにしましょう.

ここで,GPU上で実行することのできる関数オブジェクトは(ラムダ式)はrestrict制限子によってGPUで実行可能であることが修飾されていなくてはなりません.

AMPは、関数オブジェクトで渡されるカーネルに関して、restrict制限子の中でコードを書くための多くのクラスをサポートしています。

では、このparallel_for_each関数を使うために必要な知識を以下の順序で説明していきます.

1.restrict制限指定子

2.アクセラレータ

3.カーネル向けコンテナ及び関連クラス


restrict制限指定子


restrict(cpu)

通常の関数


restrict(amp)

GPUで実行可能であることを示す.

以下の条件を満たさない場合,コンパイルエラーとなる.

基本要件

1.restrict(amp)関数はrestrict(amp)関数のみ呼び出すことができる

2.int、unsigned int、float、doubleの変数と、これらの型だけを含むクラスと構造体のみを利用することができる.また,ポインタ変数もこれらの型のみをサポートする

禁止事項

1.restrict(amp)ラムダ式の参照キャプチャの利用

2.volatile,staticキーワード

3.仮想関数

4.goto文

5.dynamic_cast

6.typeid

7.可変引数関数

8.グローバル変数へのアクセス

9.再帰

まあif,forが使えるだけC++11のconstexpr関数よりはマシですね.

GPU上で実行されるため,CPUのメモリへのアクセスを伴う操作は禁止されています.

また,GPUのアーキテクチャ上シェーダで実現が困難な操作も禁止されています,


アクセラレータ

GPUデバイスの事.

利用可能なアクセラレータは以下のようにして取得できます.

#include<iostream>

#include<algorithm>
#include<vector>

using concurrency::accelerator;

std::vector<accelerator> findAccelerators(){
std::vector<accelerator> accels;
accels = accelerator::get_all();

for(int i=0; i<accels.size(); i++){
std::wcout << i+1 << "th device = " << accels[i].get_description() << "\n";
}

//emulatorのアクセラレータを削除します
accels.erase(std::remove_if(accels.begin(),accels.end(),[](accelerator& accel){return accel.get_is_emulated();}),accels.end());

return accels;
}

オンボードやエミュレータ,グラフィックボード含めてマシン全体で利用可能なGPUをget_all関数で取得できます.

エミュレータはデバッグモードでしか動かない? 正体不明のデバイスなので私はいつも除外しています.get_is_emulated関数で取得できます.

取得したデバイスの情報は以下のような関数等を用いて取得できます.

void getAccelDiscription(const accelerator& accel){

std::wcout << "accelerator: "<< accel.get_description() << std::endl;
std::cout << "version of the accelerator: " << accel.get_version() << std::endl;
std::cout << "memory: " << accel.get_dedicated_memory()/1024./1000. << " [GB]" << std::endl;;
std::cout << "is supporting double precision: " << (accel.get_supports_double_precision() ? "yes" : "no") << std::endl;
std::cout << "is attached to a display: " << (accel.get_has_display() ? "yes" : "no") << std::endl;
std::cout << "is supporting cpu shared memory: " << (accel.get_supports_cpu_shared_memory() ? "yes" : "no") << std::endl;
return;
}

倍精度浮動点小数をサポートしているか,CPUとの共有メモリをサポートしているか,等は処理において重要な情報と言えるかと.

特に共有メモリは同期のオーバーヘッドにおいて強いアドバンテージを持つと思います.

メモリが最大のアクセラレータを選ぶのはこんな感じの記述になりますね.

std::vector<accelerator>::iterator getBiggestMemoryAccelerator(std::vector<accelerator>& accels){

return std::max_element(accels.begin(),accels.end(),[](const accelerator& rhs,const accelerator& lhs){return rhs.get_dedicated_memory() < lhs.get_dedicated_memory();});
}


カーネル向けコンテナ


CPUメモリのリソースのGPUでの利用方法

parallel_for_each関数内ではメモリにアクセスできません.

渡す関数オブジェクトの中にコピーするなり,ラムダ式でコピーキャプチャを行うなりしてparallel_for_each内でGPUメモリへのコピーを行う必要があります.

parallel_for_eachは優秀ですので,デバイスをまたいだ処理でも同期処理をしているかのように振舞ってくれます.つまり,CPU側でスレッドと止めるとか,GPUとCPU間のメモリの同期を待つとか,そういったことを気にした処理を書く必要はありません.

ただし,parallel_for_each内で実行されるカーネルの順序は保証されていません.このあたりは普通の並列処理と一緒ですね.

基本的な変数はラムダ式のコピーキャプチャで事足りますが,配列の場合はポインタが使えませんので一個一個キャプチャするか,あるいは専用のAMPのコンテナを使う必要があります.

concurrency::arrayがそれに該当します.

このコンテナは普通の配列と同じように扱えますが,GPUのメモリに要素が確保されます.

template <

typename _Value_type,
int _Rank
>
friend class array;

第一テンプレート引数は型,第二引数は「次元」です.要素数ではありません.

std::arrayと混同しないように.

また,各コンテナーをGPU,CPU両方で利用可能にするラッパーコンテナであるarray_viewというものが存在します.

template <

typename _Value_type,
int _Rank = 1
>
class array_view;

このラッパーコンテナはメモリを参照し,その仲介をなすコンテナとして機能します.

array_viewに対する変更はarray_viewの同期関数によってメモリの参照元に同期されます.またデストラクタや特定の条件でも同期が行われるようです.

このラッパーコンテナを利用することで,CPUのリソースをGPUで,GPUのリソースをCPUで扱うことができます.

例えば,GPUのメモリに確保されたconcurrency::arrayに対してCPU側からアクセスするコードは以下のようなものが考えられます.

#include"amp.h"

#include<array>
#include<iostream>

template<class T,int dim,class F>
void accessArray(concurrency::array<T,dim>& vGArray,F&& function){
concurrency::array_view<T,dim> vGArrayView = vGArray; //concurrency::arrayのラッパーを作成
function(vGArrayView); //array_viewはcpu側からアクセス可能
}

template<class T, int dim>
std::unique_ptr<concurrency::array<T, dim>> createArray(const concurrency::accelerator& accel, int size) {
return std::make_unique<concurrency::array<T, dim>>(size, accel.get_default_view());
}

int main() {
constexpr int dim = 1;
const int size = 100;

concurrency::accelerator accel = *getBiggestMemoryAccelerator(findAccelerators());

auto vGArray(createArray<int, dim>(accel, size));

accessArray(*vGArray, [&](auto& _array) {
for (int i = 0; i<size; i++)_array[i] = i;
});

return 0;
}

amp関数の制限上,std::size_tはconcurrency名前空間で基本的に使われていないことも留意しておくと良いでしょう.

array_viewはコピーキャプチャを用いてGPUに転送することができます.

CPU側で確保した配列をカーネルで使いたい場合は、以下のようなコードになるでしょう.

int main() {

constexpr int dim = 1;
const int size = 100;
std::array<int,size> arr;

concurrency::accelerator accel = *getBiggestMemoryAccelerator(findAccelerators());
concurrency::extent<dim> ex;
ex[0] = size;

concurrency::array_view<int, dim> view(size,reinterpret_cast<int*>(&arr[0])); //iteratorに対応していないので仕方なくbegin()ではなくポインタをキャストしてます
parallel_for_each(accel.get_default_view(),
ex,
[=](concurrency::index<dim> gindex) restrict(amp) {
view[gindex] = 114514; //array_viewはコピーキャプチャ可能
}
);

view.synchronize(); //メモリ参照元と同期します

for (int i = 0; i<size; i++) {
std::cout << arr[i] << ",";
}
std::cout << std::endl;

return 0;
}

全ての要素を適当な数字で初期化しています.


indexクラス

template <

int _Rank
>
class index;

N次元空間の要素を一意に参照するベクタークラスです.

このクラスはparallel_for_eachに渡されるカーネルの引数になるわけですが,このインデックスがスレッド番号と結びつけられ,データの並列の処理が実現されるわけです.

多次元のarray_viewにindexを用いてアクセスするコードは以下のようになります.


template<class T,int Rank,typename... Args>
T& accessArrayByIndex(const concurrency::array_view<T,Rank>& a,Args... indexes) restrict(amp)
{
static_assert(sizeof...(indexes) == Rank,"number of index is incorrect");

concurrency::index<Rank> idx(indexes...);
return a[idx];
}

int main(void){
constexpr int COLS=6,ROWS=4;
std::array<std::array<float,COLS>,ROWS> data={
1,2,3,4,5,6,
7,8,9,10,11,12,
1,2,3,4,5,6,
7,8,9,10,11,12
};

concurrency::array_view<float,2> data_view(ROWS,COLS,reinterpret_cast<float*>(&data[0][0]));
//print data_view[3][2] by index<2>;
std::cout << accessArrayByIndex(data_view,3,2);
std::cout << std::endl;
return 0;
}

次元にかかわらず、[]オペレータで一意に要素を指定できるのが強みなわけです.

さて,この辺りで大体基本的な説明は終わりました.

このくらい分かれば並列処理はたいてい実現できると思います.

最後に, 「タイル」と呼ばれる機能を使って, 爆速で画像処理をやってみますね.


応用 タイルによる爆速画像処理

タイルとは, スレッドを四角形のサブセットに分割する機能です.

多次元でもできるかと思いますが, 二次元, 例えば画像処理で考えるとわかりやすいと思います.

例えば9*9ピクセルの画像があったとして, その画像を3*3の4つのサブセットに分割して処理をします.

このサブセット化によって何が得なのかというと, そのサブセット単位で静的変数を扱うことが出来, 同じ計算をしなくて済むようになったりします.

実際に画像を適切なサイズのサブセットに分割して, そのサブセット単位で色の平均をとる処理をしてみましょう.

すなわち, モザイク処理です.

1つのサブセットがモザイク後の画像の大きなピクセルに相当します.


index/tiled_index_modules.hpp

template<typename T,int TILE_COLS, int TILE_ROWS>

std::unique_ptr<T[]> convolutionCalculateAverage(T* data, int rows, int cols,const concurrency::accelerator& accel)
{
std::unique_ptr<T[]> average(new T[rows*cols]);

concurrency::array_view<T, 2> data_view(rows, cols, data);
concurrency::array_view<float, 2> average_view(rows, cols, reinterpret_cast<float*>(average.get()));

std::cout << "\n-------------------parallel calculation-----------------" << std::endl;
std::cout << "rows/cols " << rows << "/" << cols << std::endl;

average_view.discard_data();
parallel_for_each(
data_view.get_extent().tile<TILE_ROWS, TILE_COLS>(),
[=](concurrency::tiled_index<TILE_ROWS, TILE_COLS> idx) restrict(amp) {
tile_static T nums[TILE_ROWS][TILE_COLS];
nums[idx.local[1]][idx.local[0]] = data_view[idx.global];
idx.barrier.wait();
T sum=0;
for (int i = 0; i<TILE_ROWS; i++) {
for (int j = 0; j<TILE_COLS; j++) {
sum += nums[i][j];
}
}
average_view[idx.global] = sum / static_cast<T>(TILE_ROWS*TILE_COLS);
}
);
average_view.synchronize();

return std::move(average);
}



main.cpp

#include"opencv_include.h" //これはopencvの環境構築用のヘッダです

#include"amp.h"
#include"index/tiled_index_modules.hpp"
#include<iostream>
#include<chrono>

void image_processing_test(concurrency::accelerator& accel)
{
cv::Mat input;
cv::Mat_<float> gray,gray_cpu;
input = cv::imread("image_middle.jpg",cv::IMREAD_GRAYSCALE);
input.convertTo(gray, CV_32FC1);
input.convertTo(gray_cpu, CV_32FC1);

for (int rows = 0; rows < input.rows; rows++) {
for (int cols = 0; cols < input.cols; cols++) {
gray.at<float>(rows, cols) /= 255.;
}
}
std::cout << gray.elemSize1() << std::endl;
std::cout << gray.channels() << std::endl;
std::cout << gray.step << "/" << gray.elemSize() * gray.cols << std::endl;
std::cout << gray.isContinuous() << std::endl;

constexpr int convolution_size = 15;
{
std::chrono::time_point<std::chrono::system_clock> now = std::chrono::system_clock::now();
for (int rows = 0; rows < gray.rows; rows++) {
for (int cols = 0; cols < gray.cols; cols++) {

int sum;

for (int y = -convolution_size; y <= convolution_size; y++) {
for (int x = -convolution_size; x <= convolution_size; x++) {
if (rows + y >= 0 && rows + y < gray.rows && cols + x >= 0 && cols + x < gray.cols)
sum += gray.data[(rows + y) * gray.step + (cols + x) * gray.elemSize()];
else sum += gray.data[rows * gray.step + cols * gray.elemSize()];
}
}
gray_cpu.data[rows * gray.step + cols * gray.elemSize()] = sum / pow(2 * convolution_size + 1, 2);

}
}
std::chrono::time_point<std::chrono::system_clock> after = std::chrono::system_clock::now();
std::cout << "----------------cpu calculation succeeded---------------" << std::endl;
std::chrono::duration<double> diff = after - now;
std::cout << "score " << diff.count() << "[s]" << std::endl;
}

{
std::chrono::time_point<std::chrono::system_clock> now = std::chrono::system_clock::now();
auto result = convolutionCalculateAverage<float, convolution_size*2, convolution_size*2>(reinterpret_cast<float*>(&gray.data[0]), input.rows, input.cols, accel);
std::chrono::time_point<std::chrono::system_clock> after = std::chrono::system_clock::now();

std::cout << "----------------gpu calculation succeeded---------------" << std::endl;

std::chrono::duration<double> diff = after - now;

std::cout << "score " << diff.count() << "[s]" << std::endl;

for (int rows = 0; rows < input.rows; rows++) {
for (int cols = 0; cols < input.cols; cols++) {
gray.at<float>(rows, cols) = result[rows*input.cols + cols];
}
}
}

cv::namedWindow("window", CV_WINDOW_AUTOSIZE);
cv::namedWindow("window2", CV_WINDOW_AUTOSIZE);

cv::imshow("window", gray);
cv::imshow("window2", input);

cv::waitKey(0);
}


cpuとgpu両方で計算して, スコアを表示したりしてます.

色々説明してない関数, クラスがありますが, 疲れちゃったので開発元のmicrosoftさんに聞いてください.

ただし,cv::Matクラスにおいて,配列の次元数を2次元にするためグレースケールにしていること,

unsigned charはC++ AMPは対応していないので, 画素のサイズをfloatに変換していることに注意してください.

実行結果は以下のように.

変換前

変換後

出力

0th device = NVIDIA GeForce GTX 670

1th device = Microsoft Basic Render Driver
2th device = Software Adapter
3th device = CPU accelerator
accelerator: NVIDIA GeForce GTX 670
version of the accelerator: 720896
memory: 1.98681 [GB]
is supporting double precision: yes
is attached to a display: yes
is supporting cpu shared memory: yes
4
1
5760/5760
1
----------------cpu calculation succeeded---------------
score 1.94933[s]
-------------------parallel calculation----------------
rows/cols 810/1440
----------------gpu calculation succeeded---------------
score 0.0314382[s]

以上です.

別の環境(GTX 630M)でも試しましたが, CPUの大体50~60倍のスコアに落ち着くようです.

画像は侵略!イカ娘より引用しました.

この記事で使用したソースコードなどはgithubで公開しています.

Riyaaaaa/Gpu-Accelerated-Cpp

明日はmitsutaka-takedaさんです.よろしくおねがいします!