書籍レビュー: 「新版 独習C#」
「新版 独習C#」 山田祥寛 : 翔泳社 ISBN978-4-7981-5382-7
タイトル通り、C#の教本/自習書です。
僕は根っからのC++屋なもんで、ちょっとしたアプリ/ライブラリの類はC++で書いちまう。けどもGUIを伴うアプリあるいは.NETから使えるC++ライブラリのテストとなるとC#もそこそこ書けないと何かと不便です。
C#は15年も前のver.1.1の頃に一冊だけ買い求めました。基本的な文法/構文は現在に至るまで大きく変わることがなかったので古めかしい教本でもさほどの不便を感じてはいなかったのですが、ミョーにまだるっこい、遠回しなC#コードを書いてるんじゃないか、きょうびのC#ならもっと簡潔でエレガントに表現できるんじゃないかと思い始めていたのです。
で、モダンなC#を知ろうと言語仕様書を探してみると、公式のドキュメントはECMA-334でしてこれがC#5留まり、(2017時点で)最新のC#7にぜんぜん追いついていませんし、Microsoftの文書にしてもC#5に加えてC#6,7の差分の形でしか読むことができません。
おベンキョしよかと資料漁ってみたらばそんな有様、モヤモヤしてるところで「C#7の教本欲しくない? レビュー書いてもらうけど」のオファーが目に留まった次第です。
「新版 独習C#」ハプログラミング言語の教本としてはオーソドックスな章立てで構成されています。1～4章でイントロから始まり変数・型・演算子そして制御構文、ここまで習得すれば簡単なコンソールアプリなら作れます。
つぎにひとつメの山が訪れます。近頃の言語はどれも大抵そうなのですが、C#ってば文法/構文はシンプルにしておいて力不足な部分をライブラリで補っているように思えます。5,6章は主要かつ基本的な標準ライブラリのいくつかとコレクション(コンテナ)の解説、ココきっちり押さえておくとのちのち何かと役に立つ、ってゆーかアイデアをコードで表現するのに楽ができるはず。
7,8,9章がふたつメの山場:オブジェクト指向。Cなんかの手続き型言語に染まったビギナには考え方の大きな転換を迫られるでしょうね。手続き型だと「関数にデータを食わせて加工する」を繰り返して目的の機能を実現してたことでしょう。これがオブジェクト指向だと「データに命令して新たなデータを手に入れる」の繰り返しとでも申しましょうか、データと操作の主客が交代するよな感じになるんです。ココ躓きやすいから注意ね。
10章以降はかなり凝ったワザの数々、WPF/UWP等によるGUIアプリの実装にはほぼ必須になるから覚悟されたし。
本書の章立てを順にざっくり追いかけると以上のような雰囲気です。アタマから着々とこなしていけば、あとはライブラリ/フレームワークの使い方を学べばGUIアプリが作れるようになりますよ。
てなわけで「新版 独習C#」、とりあえず動けばいいやじゃなく基礎からきっちり学ぶべく並べられた、ストレスのない構成です...なにより日本語ですし。もちろん原本から学ぶのが間違いないのは間違いないけど...ECMA-334は当然英語、僕にはやっぱりストレスフルですわ。
本書が僕的になにより有難いのは、もくじのタイトルにC#6,7で追加された項であることが明記されていること。基本はそこそこ理解してるけど新機能には疎いんですよね。
on-memory database の検索 : SELECT JOIN ～
SELECT WHERE の応用(?)として SELECT JOIN の試作。
SELECT * FROM A, B JOIN A.id = B.id WHERE A.weight < 0.1;
ふたつの配列 A[Na], B[Nb] があって、A と B の直積すなわち A[x]. B[y] (x = 0..Na-1, y = 0..Nb-1) のすべてのペア(Na*Nb通り)に対し、A[x].id == B[y].id かつ A[x].weight < 0.1 を満たす x,y の組を列挙する、と。
直積ですから x,y に対する 二重のfor-loop で全組み合わせが作れます。CUDAではスレッド数を決定する grid, block が三次元(dim3)なので三重のloopまでならするっと書けます。
/*
SELECT * FROM arecord, brecord
JOIN arecord.id = brecord.id
WHERE arecord.weight < 0.1
*/
// nvcc fake_DB.cpp --expt-extended-lambda
// CUDA 8.0
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
struct Arecord {
int id;
float weight;
};
struct Brecord {
int id;
float height;
};
template<typename Join, typename Where>
__global__ void kernel_select(
const Arecord* arecords, unsigned int asize,
const Brecord* brecords, unsigned int bsize,
Join jpredicate, Where wpredicate,
int* count, uint2* indices) {
unsigned int sx = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int sy = blockDim.y * blockIdx.y + threadIdx.y;
if ( sx ==0 && sy == 0 ) *count = 0;
__syncthreads();
for ( unsigned int x = sx; x < asize; x += gridDim.x * blockDim.x ) {
for ( unsigned int y = sy; y < bsize; y += gridDim.y * blockDim.y ) {
if ( jpredicate(arecords[x],brecords[y]) &&
wpredicate(arecords[x],brecords[y]) ) {
indices[atomicAdd(count, 1)] = uint2{x,y};
}
}
}
}
#include <random>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
int main() {
const int Na = 100;
const int Nb = 100;
Arecord arecords[Na];
Brecord brecords[Nb];
// fill records[] with random numbers
mt19937 gen;
uniform_int_distribution<int> idist(1,10);
uniform_real_distribution<float> fdist;
generate_n(arecords, Na, [&]() { return Arecord{ idist(gen),fdist(gen) }; });
generate_n(brecords, Nb, [&]() { return Brecord{ idist(gen),fdist(gen) }; });
// result
uint2 indices[Na*Nb];
int count;
// allocate device-mem.
Arecord* dev_arecords;
cudaMalloc(&dev_arecords, Na*sizeof(Arecord));
Brecord* dev_brecords;
cudaMalloc(&dev_brecords, Nb*sizeof(Brecord));
uint2* dev_indices;
cudaMalloc(&dev_indices, Na*Nb * sizeof(uint2));
int* dev_count;
cudaMalloc(&dev_count, sizeof(int));
// copy records from HOST to DEVICE
cudaMemcpy(dev_arecords, arecords, Na*sizeof(Arecord), cudaMemcpyHostToDevice);
cudaMemcpy(dev_brecords, brecords, Nb*sizeof(Brecord), cudaMemcpyHostToDevice);
auto join = [] __device__(const Arecord& arec, const Brecord& brec) { return arec.id == brec.id; };
auto where = [] __device__(const Arecord& arec, const Brecord& brec) { return arec.weight < 0.1f; };
dim3 grid { 2, 2 };
dim3 block { 38, 8 };
kernel_select<<<grid, block>>>(
dev_arecords, Na, // FROM arecords
dev_brecords, Nb, // ,brecords
join, // JOIN arecords.id = brecords.id
where, // WHERE arecords.weight < 0.1
dev_count, dev_indices); // SELECTed int-pair
// copy result from DEVICE to HOST.
cudaMemcpy(&count, dev_count, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(indices, dev_indices, count*sizeof(uint2), cudaMemcpyDeviceToHost);
//*/
for (int i = 0; i < count; ++i) {
cout << left
<< setw(3) << indices[i].x << ":"
<< setw(10) << brecords[indices[i].y].height << " , "
<< setw(3) << indices[i].y << ":"
<< setw(10) << arecords[indices[i].x].weight << endl;
}
//*/
cout << count << " records found." << endl;
cudaFree(dev_arecords);
cudaFree(dev_brecords);
cudaFree(dev_count);
cudaFree(dev_indices);
cudaDeviceReset();
}
Device-memoryの制限がよりいっそうキツくなります。Device-memoryに載せる配列がふたつになりますし、検索結果(intのペア)の格納域がNa*Nb に比例して増えますもん。
on-memory database の検索 : SELECT WHERE ～
週末の退屈しのぎにコード書き。
struct record {
float height;
float weight;
};
なんてな struct の配列 records[N] から weight < 0.1 を満たすものを列挙します。
SQLでいうところの
SELECT * FROM records WHERE weight < 0.1;
ですわな。
これをCUDAでやるとなると records[N] および 検索結果を格納するindex列 int indices[N] はDevice上に置くことになります。
WHERE 節に与える条件式は device-lambda 使って楽しました。
// nvcc fake_DB.cpp --expt-extended-lambda
// CUDA 8.0
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
struct record {
float height;
float weight;
};
template<typename Where>
__global__ void kernel_select(
const record* records, unsigned int size,
Where condition, int* count, int* indices) {
unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
if (i == 0) *count = 0; // 一回だけ、 *count を 0-clear
__syncthreads();
while (i < size) {
// i番目のrecordが WHERE条件を満たしているなら
// indicesにiを追加
if ( condition(records[i]) ) {
indices[atomicAdd(count, 1)] = i;
}
i += gridDim.x * blockDim.x;
}
}
#include <random>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
int main() {
const int N = 100;
record records[N];
// fill records[] with random numbers
mt19937 gen;
uniform_real_distribution<float> dist;
generate_n(records, N, [&]() { return record{ dist(gen),dist(gen) }; });
// result
int indices[N];
int count;
// allocate device-mem.
record* dev_records;
cudaMalloc(&dev_records, N * sizeof(record));
int* dev_indices;
cudaMalloc(&dev_indices, N * sizeof(int));
int* dev_count;
cudaMalloc(&dev_count, sizeof(int));
// copy records from HOST to DEVICE
cudaMemcpy(dev_records, records, N * sizeof(record), cudaMemcpyHostToDevice);
// SELECT * FROM dev_records
// WHERE weight < 0.1
auto where = [] __device__(const record& rec) { return rec.weight < 0.1f; };
dim3 grid = (N + 255) / 256;
dim3 block = 256;
kernel_select <<<grid, block >> >(
dev_records, N,
where,
dev_count, dev_indices);
// copy result from DEVICE to HOST.
cudaMemcpy(&count, dev_count, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(indices, dev_indices, count * sizeof(int), cudaMemcpyDeviceToHost);
//*/
for (int i = 0; i < count; ++i) {
cout << left
<< setw(10) << records[indices[i]].height << " : "
<< setw(10) << records[indices[i]].weight << endl;
}
//*/
cout << count << " records found." << endl;
cudaFree(dev_records);
cudaFree(dev_count);
cudaFree(dev_indices);
cudaDeviceReset();
}
メモリの乏しいDeviceでは扱えるrecord数をそんなに大きくはできません。まぢに使うとなればrecords[N]を分割し、Deviceに乗せては検索を繰り返すことになるでしょうね。
Grid-stride loop
NVIDIAの開発系BlogになるほどなーなTips: Grid-stride loop を見つけたので軽く解説。 (元ネタはコチラ)
毎度おなじみ SAXPY : ベクトルの積和演算 Y[i] = a * X[i] + Y[i] (i = 0,1 ... n-1) これをCPUでフツーに行うとき、iについてくるくる回すloopで実装しますわな。
void saxpy(int n, float a, const float* x, float* y) {
for ( int i = 0; i < n; ++i ) {
y[i] = a * x[i] + y[i];
}
}
かたやこいつをGPU(CUDA)でやるときの常套手段は "iによるn回のloop" を "n個のスレッド" に置き換えることで コードからloopを取っ払います。
__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
int i = blockDim.x * blockIdx.x + threadIdx.x;
if ( i < n ) {
y[i] = a * x[i] + y[i];
}
}
void saxpy(int n, float a, const float* x, float* y) {
kernel_saxpy<<<(n+255)/256, 256>>>(n, a, x, y);
}
Grid-stride loopは両者の中間みたいなやりかたで、for-loopを温存します。
__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
for ( int i = blockDim.x * blockIdx.x + threadIdx.x;
i < n;
i += gridDim.x * blockDim.x ) {
y[i] = a * x[i] + y[i];
}
}
このとき iの初期値はスレッドの通し番号、loopごとの増分はスレッドの総数すなわち gridDim.x * blockDim.x です。するってーとスレッド総数が Tのとき、
スレッド#0 : i = 0, T , 2T ... スレッド#1 : i = 1, T+1, 2T+1 ... スレッド#2 : i = 2, T+2, 2T+2 ... (以下同文)
に対して y[i] を計算してくれる、と。
こうしておくとなにがオイシイかというと、kernel-call時のスレッド数を自由に指定できます。grid * blockがnより小さくてもちゃんと長さnのベクトルについて計算してくれる。
なのでたとえばGPUが持ってるプロセッサ数に応じてgrid/block数をチューニングできます。
void saxpy(int n, float a, const float* x, float* y) {
int numSMs;
cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
// SMあたり32*256=8192個のスレッドを起こして計算する
kernel_saxpy<<<numSMs*32, 256>>>(n, a, x, y);
}
極端な話 kernel_saxpy<<<1, 1>>>(n, a, x, y); でも問題ありません。このときスレッド一個で処理しますから多数のスレッドが絡む問題を起こさないので本来得られるべき結果が出てくるわけで、デバッグ時の検証に重宝しますね。
加えて kernel内で'prinf'したとき出力がお行儀よく並びますからデバッグが楽ですよね。
おためしに書いてたコードがこちらになります:
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
for ( int i = blockDim.x * blockIdx.x + threadIdx.x;
i < n;
i += gridDim.x * blockDim.x ) {
y[i] = a * x[i] + y[i];
}
}
#include <iostream>
#include <random>
#include <algorithm>
#include <ctime>
int main() {
using namespace std;
const int N = 100000;
float* x;
float* y0;
float* y;
float a;
cudaMallocManaged(&x, N*sizeof(float));
cudaMallocManaged(&y, N*sizeof(float));
cudaMallocManaged(&y0, N*sizeof(float));
mt19937 gen(time(nullptr));
uniform_real_distribution<float> dist;
auto rand = [&]() { return dist(gen); };
generate_n(x, N, rand);
generate_n(y, N, rand);
copy(y, y+N, y0);
a = rand();
// 単一スレッドで計算する
kernel_saxpy<<<1,1>>>(N, a, x, y0);
int numSMs;
cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
cout << numSMs << " SMs" << endl;
// SMあたり32*256=8192個のスレッドを起こして計算する
kernel_saxpy<<<numSMs*32, 256>>>(N, a, x, y);
cudaDeviceSynchronize();
if ( equal(y0, y0+N, y) ) {
cout << "ok." << endl;
}
cudaFree(x);
cudaFree(y);
cudaFree(y0);
cudaDeviceReset();
}
転置しながらコピーする
cuBLAS が扱う行列要するに二次元配列のメモリレイアウトはFortranとC/C++では大きく異なります。 行列
│ 1 2 3 │ │ 4 5 6 │
がメモリ上にリニアに展開されたとき、
Fortranでは 1 4 2 5 3 6 の順、上から下/左から右に並ぶのに対し、
C/C++ では 1 2 3 4 5 6 の順、左から右/上から下に並びます。
上から下(行:row方向)が左から右(列:column方向)より優先するFortran流レイアウトをrow-major、反対のC/C++流レイアウトをcolumn-majorっていいます。
cuBLASの基となった本家BLASはもともとFortranで実装されたライブラリであるがため、行列表現はFortranに倣ってrow-majorになってます。したがって上記の行列をC/C++の二次元配列で表現すると、
// 列 行
float A[3][2] = {
{ 1.0f, 4.0f },
{ 2.0f, 5.0f },
{ 3.0f, 6.0f }
};
となります。行列を転置(行/列を交換)せにゃならんのです。
これがどうにも居心地よろしくありません。Host上ではcolumn-major、cuBLASが扱うDevice上ではrow-majorで表現し、Host-Device間のコピーついでに転置する transpose_copy をこしらえてみました。
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <iostream>
template<typename T>
cudaError_t transpose_copy( T* dst, int ldDst, // コピー元
const T* src, int ldSrc, // コピー元
int m, int n, // 行, 列
cudaStream_t stream = nullptr) {
cudaError_t status = cudaSuccess;
for ( int i = 0; i < n && status == cudaSuccess; ++i ) {
status = cudaMemcpy2DAsync(dst, sizeof(T),
src, sizeof(T)*ldSrc,
sizeof(T), m,
cudaMemcpyDefault, stream);
++src;
dst += ldDst;
}
return status;
}
int main() {
const int LD = 8;
// 2行3列のA と
float A[LD][LD] = {
{ 1.0f, 1.1f, 1.2f },
{ 1.3f, 1.4f, 1.5f }
};
// 3行2列のB の積を
float B[LD][LD] = {
{ 2.0f, 2.1f },
{ 2.2f, 2.3f },
{ 2.4f, 2.5f }
};
// 2行2列のC に求める
float C[LD][LD];
float* dA;
float* dB;
float* dC;
cudaMalloc(&dA, LD*LD*sizeof(float));
cudaMalloc(&dB, LD*LD*sizeof(float));
cudaMalloc(&dC, LD*LD*sizeof(float));
// Host->Device
transpose_copy(dA, LD, (float*)A, LD, 2, 3);
transpose_copy(dB, LD, (float*)B, LD, 3, 2);
// Cは0で埋めておく
cudaMemset(dC, 0, LD*LD*sizeof(float));
cublasHandle_t handle;
cublasCreate(&handle);
// GEMM: C = αAB + βC (α=1, β=0 なら C = AB)
float alpha = 1.0f;
float beta = 0.0f;
cublasSgemm(handle,
CUBLAS_OP_N, CUBLAS_OP_N,
2, 2, 3,
&alpha, dA, LD, dB, LD,
&beta, dC, LD );
// Device->Host
transpose_copy((float*)C, LD, dC, LD, 2, 2);
cudaDeviceSynchronize();
for ( int y = 0; y < 2; ++y ) {
for ( int x = 0; x < 2; ++x ) {
std::cout << C[y][x] << ' ';
}
std::cout << std::endl;
}
cublasDestroy(handle);
cudaFree(dA);
cudaFree(dB);
cudaFree(dC);
}
nvgraphPagerankでPagerankを求める
CodeZineにPagerankアルゴリズムの解説記事を書きました。
推移確率行列:H に対し
λI = HI
を満たす固有ベクトル(eigenvector):I がすなわち Pagerank となります。
※ Iは確率ベクトルなので総和は常に1、なので固有値(eigenvalue):λ は1.0 です。
この問題を解いてくれるのが nvGRAPH の nvgraphPagerank、CodeZineで使ったサンプルデータで nvgraphPagerank を使ってみましたよ。
※ codeの半分は COO-format から CSC-format への変換でできてます。
#include <cuda_runtime.h>
#include <nvgraph.h>
#include <iostream>
#include <algorithm>
using namespace std;
void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float*bookmark, float* probs, int iter);
template<typename T>
void dump(T* data, int n) {
T* items = new T[n];
cudaMemcpy(items, data, n*sizeof(T), cudaMemcpyDefault);
for_each( items, items+n, [](T item) { cerr << item << " "; });
cerr << endl;
delete[] items;
}
int main() {
const float r0 = 0.0f;
const float r1 = 1.0f;
const float r2 = 1.0f / 2.0f;
const float r3 = 1.0f / 3.0f;
const float r8 = 1.0f / 8.0f;
// 元ネタ: COO-format
const int nvertices = 8; // ページ数
const int nedges = 17; // リンク数
int source_indices[nedges] = // リンク元
{ 0, 0, 1, 2, 2, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7 };
int destination_indices[nedges] = // リンク先
{ 1, 2, 3, 1, 4, 1, 4, 5, 5, 6, 7, 7, 0, 4, 7, 5, 6 };
float weights[nedges] = // 各リンクの重み
{ r2, r2, r1, r2, r2, r3, r3, r3, r3, r3, r3, r1, r3, r3, r3, r2, r2 };
// リンク先が少なくともひとつあるページは 1.0f、
// リンク先のない行き止まりページは 0.0f
float bookmark[nedges];
fill_n(bookmark, nedges, 1.0f);
for ( int dest : destination_indices ) {
bookmark[dest] = 0.0f;
}
// 上記の元ネタを coo_topology に転写
nvgraphCOOTopology32I_st coo_topology;
float* coo_weights;
cudaMalloc(&coo_topology.source_indices , nedges*sizeof(int));
cudaMalloc(&coo_topology.destination_indices, nedges*sizeof(int));
cudaMalloc(&coo_weights, nedges*sizeof(float));
coo_topology.nvertices = nvertices;
coo_topology.nedges = nedges;
cudaMemcpy(coo_topology.source_indices , source_indices , nedges*sizeof(int), cudaMemcpyDefault);
cudaMemcpy(coo_topology.destination_indices, destination_indices, nedges*sizeof(int), cudaMemcpyDefault);
coo_topology.tag = NVGRAPH_UNSORTED;
cudaMemcpy(coo_weights, weights, nedges*sizeof(float), cudaMemcpyDefault);
// nvGRAPH の出番はここから
nvgraphHandle_t handle;
nvgraphCreate(&handle);
cudaDataType_t data_type = CUDA_R_32F; // 32bit-real (=float)
{
// Pagerankを求めるには CSR-formatに変換せにゃならん
nvgraphCSCTopology32I_st csc_topology;
float* csc_weights;
cudaMalloc(&csc_topology.source_indices, nedges*sizeof(int));
cudaMalloc(&csc_topology.destination_offsets, (nvertices+1)*sizeof(int));
cudaMalloc(&csc_weights, nedges*sizeof(float));
nvgraphConvertTopology(handle,
NVGRAPH_COO_32, &coo_topology, coo_weights, // in type, topology, weigts
&data_type, // type of weights
NVGRAPH_CSC_32, &csc_topology, csc_weights); // out type, topology, weigts
cerr << "CSC : " <<endl;
cerr << "nvertices = " << csc_topology.nvertices << endl;
cerr << "nedges = " << csc_topology.nedges << endl;
cerr << "source_indices = "; dump(csc_topology.source_indices , csc_topology.nedges);
cerr << "destination_offsets = "; dump(csc_topology.destination_offsets, csc_topology.nvertices+1);
cerr << "weights = "; dump(csc_weights, csc_topology.nedges);
// Pagerankを求めて出力
float probs[nvertices] = { r1, r0, r0, r0, r0, r0, r0, r0 };
pagerank(handle, &csc_topology, csc_weights, bookmark, probs, 100);
cout << "\nPagerank: \n";
for ( int i = 0; i < nvertices; ++i ) {
cout << "page[" << i << "] : " << probs[i] << endl;
}
cudaFree(csc_topology.source_indices);
cudaFree(csc_topology.destination_offsets);
cudaFree(csc_weights);
}
cudaFree(coo_topology.source_indices);
cudaFree(coo_topology.destination_indices);
cudaFree(coo_weights);
nvgraphDestroy(handle);
cudaDeviceReset();
}
void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float* bookmark, float* probs, int iter) {
nvgraphGraphDescr_t descr;
// setup
nvgraphCreateGraphDescr(handle, &descr);
// topology
nvgraphSetGraphStructure(handle, descr, topology, NVGRAPH_CSC_32);
cudaDataType_t dimT[2] = { CUDA_R_32F, CUDA_R_32F };
nvgraphAllocateVertexData(handle, descr, 2, dimT);
nvgraphAllocateEdgeData(handle, descr, 1, dimT);
// weight
nvgraphSetEdgeData(handle, descr, weights, 0);
// bookmark
nvgraphSetVertexData(handle, descr, bookmark, 0);
// initial guess
nvgraphSetVertexData(handle, descr, probs, 1);
// pagerank!!
float damping_factor = 0.85f; // 85%はリンクを辿り,残る15%は気ままにジャンプ
nvgraphPagerank(handle, descr,
0, // weight : EdgeData[0]
&damping_factor,
0, // bookmark : VertexData[0]
1, // has guess
1, // pagerank : VertexData[1]
1.0e-6f, // tolerance 誤差 : これより小さくなったら計算終了
iter // 最大繰り返し回数
);
// get result
nvgraphGetVertexData(handle, descr, probs, 1);
// teardown
nvgraphDestroyGraphDescr(handle, descr);
}
thrust::for_eachによるヒストグラム
Thrust、個人的にはかなり気に入ってるライブラリなんだけど、なんだか影が薄いんですよね。 cuFFTやcuDNNなどなど、僕らにはとても書けそうにないガチにヘビィな連中とは異なり、Thrustが提供するのは copy/search/sortといった基本的なアルゴリズム群だからでしょうか。
CUDAが、というかnvccが device lambda をサポートしてくれたおかげで、Thrustの使い心地がぐっと良くなりました。その一例。
ヒストグラム(度数分布表)を作ります。たくさんの数値データがあって、それらのとる値の範囲を等間隔に刻んだスロット列を用意します。各数値をスロットに振り分けたとき、各スロットにいくつずつ納まるかってやつ。
アルゴリズムはこんなカンジ:
float score[N]; // 数値データ列(input)
int hist[S]; // ヒストグラム(output) 初期値はall-0
for ( score 内の各数値 val について ) {
if ( valがスロットiの範囲にあるなら } {
++hist[i};
}
}
これをkernelで実装するのはちっとも難しくありません。けどもたったこれだけのためにわざわざkernel書くのが面倒なわけで、Thrustでさらっと実装してみました。
#include <thrust/device_vector.h>
#include <thrust/for_each.h>
#include <curand.h>
#include <iterator>
#include <iostream>
#include <ctime>
int main() {
using namespace std;
const int N = 100000;
thrust::device_vector<float> dscore(N);
thrust::device_vector<int> dhist(101, 0);
{ // 平均 50.0, 標準偏差 15.0 の正規分布乱数で dscoreを埋める
curandGenerator_t gen;
curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT);
curandSetPseudoRandomGeneratorSeed(gen, static_cast<unsigned long long>(time(nullptr)));
curandGenerateNormal(gen, dscore.data().get(), N, 50.30f, 15.0f);
curandDestroyGenerator(gen);
}
// dscoreの各値を四捨五入し、0以上/100以下ならばヒストグラムに加える
int* histPtr = dhist.data().get(); // for_each内のlambdaがコレをキャプチャする
thrust::for_each(begin(dscore), end(dscore),
[=] __device__ (float val) -> void {
int i = static_cast<int>(val+0.5f); // 四捨五入して
if ( i >= 0 && i <= 100 ) { // 0以上100以下なら、当該スロットを+1する
atomicAdd(histPtr+i, 1); // ++histPtr[i] だと data-race起こすよきっと!
}
});
//*/
for ( int item : dhist ) {
cout << item << endl;
}
//*/
}
thrust::for_each(こっから, ここまで, なんかする) の"なんかする"を devie lambda で与えています。ヒストグラムを求める領域:dhistのポインタをlambda式の値キャプチャで食わせています。
※ device lambdaでキャプチャできるのは値キャプチャだけ。参照キャプチャするとコンパイル・エラーとなります。Hostにある変数の参照(つまりポインタ)をDevice側に持ち込んでも意味ないもんね。
結果をCSVにリダイレクトし、excelに食わせたものがこちらになります。
