というわけで続きである。今回はIsingモデルを書く。openGLとかで動画を眺めたりはしない。
半年ほどまえ、仲間内でIsingモデルをどの程度早く走らせることができるか競ったことがあり、そのときはOpenMPとかを使って6スレッドくらいで並列化して3072 * 3072のスピン全体の反転100回で7.5秒というのが記録だった。今回はCUDAで同様のことをしてみよう。
当時のコードを見ると少しｳｯとなるところがあることを思うと、多分数カ月か半年くらいしてこのエントリを見たらｳｯとなるのだろう。それくらいのノリでやっていく。
Ising model
なんの変哲もない2次元Isingモデルをメトロポリスモンテカルロ法でサンプリングする。ハミルトニアンは
で、(i,j)はノイマン近傍とする。
エネルギーの計算
経験的にexpの計算はそこそこ重い。しかし今回は、単純なIsingモデルであることが幸いしてexpが必要な計算は先にやってしまっておくことができる。隣り合うスピンの位置はエネルギーとは無関係なので、実質的に隣接するスピンのパターンは、自身と同じ向きのスピンが0,1,2,3,4個あるという5通りに絞ることができる。全体に磁場をかけることを考慮すると変更する対象のスピンの上下を考える必要があるが、それでもたったの10通りだ。しかも、次の候補として提案される状態は常にそのときのスピンが反転した状態しかない。なので、計算しないといけないパターンは10通りにしかならないのである。
というわけで、計算を始める前にメトロポリスモンテカルロ法のaccept確率に相当する
を決めてしまおう。
const float exp_dE[10] = { // case {neighbors}, center std::exp(beta * ( 4*J + 2*H)), // {up, up, up, up}, down std::exp(beta * ( 4*J - 2*H)), // {dn, dn, dn, dn}, up std::exp(beta * ( 2*J + 2*H)), // {up, up, up, dn}, down std::exp(beta * ( 2*J - 2*H)), // {dn, dn, dn, up}, up std::exp(beta * ( 0*J + 2*H)), // {up, up, dn, dn}, down std::exp(beta * ( 0*J - 2*H)), // {dn, dn, up, up}, up std::exp(beta * (-2*J + 2*H)), // {up, dn, dn, dn}, down std::exp(beta * (-2*J - 2*H)), // {dn, up, up, up}, up std::exp(beta * (-4*J + 2*H)), // {dn, dn, dn, dn}, down std::exp(beta * (-4*J - 2*H)) // {up, up, up, up}, up };
さらに、__constant__メモリはこのような用途に最適である。これはありとあらゆるthreadが必要とする情報で、しかもGPU側はスピンの反転のみを受け持つのでこの内容を変更しない。
modernなCUDAの書き方をまだ知らないので、グローバルに__constant__メモリを定義して、そこに上記の配列をコピーした。
// pre-calculated exp(dE / (kB * T)) __constant__ float exp_dE_beta[10]; int main() { const cudaError_t err_dE = cudaMemcpyToSymbol(exp_dE_beta, exp_dE, sizeof(float) * 10); assert(err_dE == 0); }
シミュレーションの途中で温度を上下などしようと思った時は、CPU側からは__constant__メモリは書き換えられるので計算しなおしたあとでまた書き込めばよいのだろう。
あといい加減cudaError_tを扱う関数なりなんなりを書くべきだと思う。
乱数
乱数の生成も非常にしんどい処理である。しかも、各threadが必要となったタイミングで乱数生成をしていると再現性が光の彼方に消え失せてしまう。そのあたりを正しく考慮して自分で乱数を発生させるのはやめたほうがいいだろう。しかし、CPUで乱数を生成して毎回cudaMemcpyしていては速度が出なかろう。
なのでcuRANDの出番だと考えた。cuRANDによってGPUでデバイス側で乱数を生成し、そのままデバイスでその乱数を使えばメモリ転送のオーバーヘッドはないだろうと予想し、cuda toolkit documentationで使い方を調べた。
curand.hをインクルードし、まず乱数を保管する配列をデバイス側に確保する。で、乱数生成器(curandGenrator_t)を作り、シードを設定する。
#include <curand.h> int main() { // その他諸々の前処理 float *random; cudaError_t err_random = cudaMalloc((void**)&random, sizeof(float) * w * h); assert(err_random == 0); // prepair cuRAND generators curandGenerator_t rng; curandStatus_t st_gen = curandCreateGenerator(&rng, CURAND_RNG_PSEUDO_DEFAULT); assert(st_gen == CURAND_STATUS_SUCCESS); curandStatus_t st_seed = curandSetPseudoRandomGeneratorSeed(rng, seed); assert(st_seed == CURAND_STATUS_SUCCESS); // 時間発展とか }
実際に乱数を生成するときは、host側から
// generate random numbers
curandStatus_t st_genrnd = curandGenerateUniform(rng, random, w * h);
assert(st_genrnd == CURAND_STATUS_SUCCESS);
のようにすればよい。これで、0.0から1.0までのfloatの乱数が生成できる。多分curandGenerateUniform以外の関数なら別の分布を生成できるはずだ。あんまり調べてないけど。
最後に
curandStatus_t destroy = curandDestroyGenerator(rng); assert(destroy == CURAND_STATUS_SUCCESS);
のようにして生成器を開放しておくのを忘れない。
実装
準備ができたところで実際に時間発展させるところを記述する。Isingモデルは上下左右のペアとしかスピンを比較しない(ノイマン型近傍)ので、市松模様の白い方のみ、黒い方のみ、という風に更新すれば、別のthreadが比較対象のスピンを知らないうちに反転させてしまうということを心配せずに好き勝手な順序で更新できる。
最初は必要な量のthreadだけを立ち上げて、更新するIndexをうまく設定しようと考えていたのだが、いかんせんそこは初心者、バグを起こしてとりあえず一番素直に実装することにした。つまり、更新してはいけないスピンを担当することになったthreadを即時撤退させるのである。
__global__ void update_field(bool* spins, const float* random, const std::size_t x_size, const std::size_t y_size, bool turn) { const int x = threadIdx.x + blockIdx.x * blockDim.x; const int y = threadIdx.y + blockIdx.y * blockDim.y; if(turn) { if((x+y)%2 == 1) return; } else { if((x+y)%2 == 0) return; } const std::size_t xdim = blockDim.x * gridDim.x; const std::size_t offset = x + y * xdim; if(offset >= x_size * y_size) return; // 周期境界条件を考慮して隣のIndexを決める const std::size_t n_offset = (y+1 < y_size) ? x + (y+1) * xdim : x; const std::size_t e_offset = (x+1 < x_size) ? (x+1) + y * xdim : y * xdim; const std::size_t s_offset = (y-1 >= 0) ? x + (y-1) * xdim : x + (y_size-1) * xdim; const std::size_t w_offset = (x-1 >= 0) ? (x-1) + y * xdim : x_size - 1 + y * xdim; const bool c = spins[offset]; // center const bool n = spins[n_offset]; // north const bool e = spins[e_offset]; // east const bool s = spins[s_offset]; // south const bool w = spins[w_offset]; // west // スピンの向きを比較し、事前に計算してあったexp(dE)にアクセスするIndexを決める std::size_t dJ = 0; if(c == n) ++dJ; if(c == e) ++dJ; if(c == s) ++dJ; if(c == w) ++dJ; const std::size_t dH = c ? 1 : 0; // 遷移確率と乱数を比較してacceptするか決める if(exp_dE_beta[dH + dJ * 2] > random[offset]) spins[offset] = (!c); return; }
メインループでは、乱数を生成して、これをturn = true, falseの2回呼び出すことで一回の時間発展とする。
その他
他にファイルからいくつかのパラメータを読み込めるようにしたり画像を出力できるようにしたが、本筋と逸れるので書かない。githubレポジトリには載っているのでどうしても見たければそちらを参照してほしい。
速度
さて、CUDAで実行しているのにCPU + openMPよりも遅かったりすると話にならない。時間を計測しよう。少しでも早く走るかのように見せるために出力はしない。実際、CPUで7秒とかの記録も出力はしていない（最適化ですべてが消え失せたらどうするんだ、という心配はあるが、そこまでコンパイラは賢くないと信じている）。
計測方法は一番脳みそを使わなくていいtimeコマンドである。「雑」の一言に尽きる。
$ time ./a.out input.dat 0.27user 0.64system 0:00.92elapsed 99%CPU (0avgtext+0avgdata 90160maxresident)k 0inputs+8outputs (0major+3829minor)pagefaults 0swaps
マジか。正直に言って殆どチューニングらしいチューニングはしていないのに（かろうじて遷移確率を先に計算しておいたのはチューニングといってよさそうだが）、1秒を切ってしまった。CPUで頑張って10秒切ったりしてたのは何だったんだ。
もちろんこの速度はblockとthreadの数に依る。thread数を何も考えずに1にしていたら12秒くらいかかった。上の結果は(16, 16)threadでの結果である。もう少しいろいろ試してみるか勉強することでもう少しは早くできる気がするが、それよりも必要な分のthreadだけを起動するとか、メモリアクセスを考えてみるといったことの方が速度が上がりそうな気もする。ちなみにこれはblockとthreadの最適な数を探索したくないことの言い訳である。
出力
流石に出力なしだと単に電力の無駄遣いなので、ちゃんと出力してみる。png画像を作ると1分くらいかかったので（多分これももっと高速化できるが面倒だった）、文字列として出力すればもう少し早くなるだろう、ということで
std::ofstream ofs("ising_traj.dat"); char *traj = new char[w * h]; for(/*main loop*/) { cudaError_t ercpy = cudaMemcpy( snapshot, spins, sizeof(bool) * w * h, cudaMemcpyDeviceToHost); assert(ercpy == 0); for(std::size_t i=0; i<w*h; ++i) traj[i] = static_cast<char>(snapshot[i]) + 48; ofs << traj << std::endl; }
のようにした。すると7秒程度になった。何かしら解析するときは画像を読むよりこっちの方が楽かもしれないので、まあこんなところにしておこう。
動画はこんな感じになる。ちなみにこれのサイズは1024 * 1024である。
今回のコードは
GitHub - ToruNiina/cuda_practice
のisingディレクトリにある。
