2009年6月22日月曜日

C++: インテル スレッディング・ビルディング・ブロックを使って簡単に並列化

このブログ記事をはてなブックマークに追加

最近はマルチコアCPUが当たり前になってきて、それを使って簡単に並列処理プログラミングができないか、頭を悩ませていたのだが、インテルのスレッディング・ビルディング・ブロック(Threading Building Blocks, TBB)が非常に良くできた技術であることを知った。

これまではMPIやら、OpenMPやら、CUDAやらをちょろちょろ手を出しつつも、どれも正直面倒だった。趣味のプログラミングではそれでも楽しいからいいのだけど、単にある処理を並列化で高速にしたいだけの場合、並列処理に関わる面倒な手続きなどは極力省きたい。また、どうしても泥臭いやり方になることが多い。例えば、Cでは簡単に書けるのだけど、C++の機能を使った処理がやりづらかったりね。

そこで、TBBが登場する。マルチコアやメニーコアプラットフォーム上での利用となるが、最近では一般のノートPCでさえマルチコアCPUが使われていることを考えると、世にある様々な並列処理環境のうち、マルチコアCPUでの並列化がもっとも望まれているだろうと思う。実際、自分もマルチコアCPU上での並列処理が簡単にできないものかとよく考えていた。それに、従来の泥臭さの残る並列化に比べて、C++テンプレートによるコーディングスタイルはスマートだ。

さて、TBBを利用するにはTBB公式サイトからダウンロードしてくればよい。オープンソース版が用意されているので、それを使って自由にプログラムを作成できる。また、このサイトにはチュートリアルやリファレンスも用意されている。さらに、日本語版のドキュメントもXLsoftのTBBサイトから利用できる。ここのチュートリアルを一通り読めば大抵の場合で利用できるようになるだろう。

今回は、「はじめてのCUDAプログラミングで分子動力学計算」で作成したC++版のなんちゃって分子動力学(MD)計算プログラムをTBBを利用して並列処理を行ってみた。非常に簡単だった。相互作用計算の関数を関数オブジェクトに変更しただけだ。TBBについて知っていれば5分もかからないかもしれない。こんなに簡単に並列化ができるのであれば、もっと早くに知りたかった。まあ、百聞は一見に如かず、今回使ったコードを下に示す。青(水色)で示した部分が書き加えたコード、赤(ピンク)で示した部分が削ったコード、水色とピンクのコードは共通のコードである。実際に修正した部分は青と赤で示したコードだが、それが全体に対して如何に少ないか分かってもらえるだろうか。

MDプログラムの説明は上述のエントリに書かれている。また、利用したTBBについては、まず、main関数でtask_scheduler_init init;として初期化を行っている。そして、相互作用の計算を行うSimpleMD::Calcメンバ関数をSimpleMDCalc関数オブジェクトとして実装し、SimpleMD::Runメンバ関数内でparallel_for(blocked_range(0, num_atoms, 500), SimpleMDCalc(num_atoms, crd, f, ene, chg));として呼び出している。blocked_rangeの第3引数で粒度(反復回数)を指定している。TBBについての詳しい説明については、XLsoftのチュートリアルを参照して欲しい。

最後に実際の計算時間を測定したので、それを示す。チャートではコアに対する実行速度向上の割合を示している。4コアのXeonが8個搭載されたSMP環境で、19,393個の疑似原子を中心から50Å以内に発生させて5ステップの計算を行った。因みにコア数の指定はtask_scheduler_init init;で行う。例えば8コアを指定する場合は、init(8);とすればよい。

single 78.18秒 2 cores 39.28秒 1.985倍 4 cores 22.67秒 3.449倍 8 cores 11.70秒 6.682倍 16 cores 7.27秒 10.784倍

TBB

さすがに複数コアを利用するとある程度はサチってしまうが、それでもこれだけ簡単に利用できることを考えると十分な性能が発揮されていると思う。

simple_md_tbb.cpp

///////////////////////////////////////////////////////////////////// // Simple Molecular Dynamics by nox, 2009.06.22 // // Perform molecular dynamics for pseudo-atoms without internal // interactions, such as bonds, angles, dihedrals. Only vdW and // Coulomb interactions. #include <iostream> #include <iomanip> #include <fstream> #include <sstream> #include <cmath> #include "tbb/task_scheduler_init.h" #include "tbb/parallel_for.h" #include "tbb/blocked_range.h" using namespace std; using namespace tbb; class SimpleMD { private: int num_steps; int num_atoms; double cap_range; double* crd; double* f; double* ene; double* chg; public: SimpleMD(int, char*); ~SimpleMD(); void ReadCrd(char*); void Center(); void Calc(); void Move(); void PrintCrd(); void Run(); double Energy() const { double sum = 0.0f; for (int i = 0; i < num_atoms; i++) sum += ene[i]; return sum / 2.0f; } double DistanceXYZ(int i, int j, int k) const { return crd[i * 3 + k] - crd[j * 3 + k]; } double Distance(int i, int j) const { double dx = crd[i * 3] - crd[j * 3]; double dy = crd[i * 3 + 1] - crd[j * 3 + 1]; double dz = crd[i * 3 + 2] - crd[j * 3 + 2]; return sqrt(dx * dx + dy * dy + dz * dz); } }; // Calculate energies and forces // Total energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i])^6] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz class SimpleMDCalc { private: int num_atoms; double* crd; double* f; double* ene; double* chg; public: SimpleMDCalc(int num_atoms_, double* crd_, double* f_, double* ene_, double* chg_) : num_atoms(num_atoms_), crd(crd_), f(f_), ene(ene_), chg(chg_) { } double DistanceXYZ(int i, int j, int k) const { return crd[i * 3 + k] - crd[j * 3 + k]; } double Distance(int i, int j) const { double dx = crd[i * 3] - crd[j * 3]; double dy = crd[i * 3 + 1] - crd[j * 3 + 1]; double dz = crd[i * 3 + 2] - crd[j * 3 + 2]; return sqrt(dx * dx + dy * dy + dz * dz); } void operator()(const blocked_range<int>& br) const { const double eps = 0.2f; const double Rij = 2.5f; double r, r_xyz, df, r12, r6, q; for (int i = br.begin(); i != br.end(); i++) { ene[i] = 0.0f; for (int j = i + 1; j < num_atoms; j++) { r = Distance(i, j); q = chg[i] * chg[j]; r12 = pow(Rij / r, 12); r6 = pow(Rij / r, 6); ene[i] += 2.0 * (eps * (r12 - 2.0f * r6) + q / r); for (int k = 0; k < 3; k++) { r_xyz = DistanceXYZ(i, j, k); df = (12 * eps / Rij) * (r12 * Rij / r - r6 * Rij / r) * r_xyz / r; df += q / pow(r, 3) * r_xyz; f[i * 3 + k] += df; f[j * 3 + k] -= df; } } } } }; SimpleMD::SimpleMD(int n, char* fname) : num_steps(n) { ReadCrd(fname); f = new double[num_atoms * 3]; ene = new double[num_atoms]; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) f[i * 3 + j] = 0.0f; ene[i] = 0.0f; } Center(); } SimpleMD::~SimpleMD() { delete[] ene; delete[] f; delete[] chg; delete[] crd; } void SimpleMD::ReadCrd(char* fname) { fstream fs(fname, ios_base::in); string line; stringstream ss; if (!fs.is_open()) { cerr << "File open error: " << fname << endl; exit(1); } getline(fs, line); cout << line << endl; getline(fs, line); ss.str(line); ss >> num_atoms; ss.clear(); getline(fs, line); ss.str(line); ss >> cap_range; ss.clear(); crd = new double[num_atoms * 3]; chg = new double[num_atoms]; for (int i = 0; i < num_atoms; i++) { getline(fs, line); ss.str(line); ss >> crd[i * 3] >> crd[i * 3 + 1] >> crd[i * 3 + 2] >> chg[i]; ss.clear(); ss.str(""); } fs.close(); } void SimpleMD::Center() { double d[3] = { 0.0f, 0.0f, 0.0f }; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) d[j] += crd[i * 3 + j]; } for (int j = 0; j < 3; j++) d[j] /= num_atoms; for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) crd[i * 3 + j] -= d[j]; } void SimpleMD::Calc() { const double eps = 0.2f; const double Rij = 2.5f; double r, r_xyz, df, r12, r6, q; for (int i = 0; i < num_atoms; i++) { ene[i] = 0.0f; for (int j = i + 1; j < num_atoms; j++) { r = Distance(i, j); q = chg[i] * chg[j]; r12 = pow(Rij / r, 12); r6 = pow(Rij / r, 6); ene[i] += 2.0 * (eps * (r12 - 2.0f * r6) + q / r); for (int k = 0; k < 3; k++) { r_xyz = DistanceXYZ(i, j, k); df = (12 * eps / Rij) * (r12 * Rij / r - r6 * Rij / r) * r_xyz / r; df += q / pow(r, 3) * r_xyz; f[i * 3 + k] += df; f[j * 3 + k] -= df; } } } } void SimpleMD::Move() { double r, dr; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) { crd[i * 3 + j] += f[i * 3 + j] * 0.01f; f[i * 3 + j] = 0.0f; } r = 0.0f; for (int j = 0; j < 3; j++) r += crd[i * 3 + j] * crd[i * 3 + j]; dr = r - cap_range * cap_range; if (dr > 0.0f) { for (int j = 0; j < 3; j++) f[i * 3 + j] -= crd[i * 3 + j] / cap_range * dr * 0.01f; } } } void SimpleMD::PrintCrd() { cout << endl; cout << num_atoms << endl; cout << cap_range << endl; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) cout << " " << fixed << setw(10) << setprecision(6) << crd[i * 3 + j]; cout << " " << fixed << setw(8) << setprecision(4) << chg[i]; cout << endl; } } void SimpleMD::Run() { for (int i = 0; i < num_steps; i++) { Calc(); parallel_for(blocked_range<int>(0, num_atoms, 500), SimpleMDCalc(num_atoms, crd, f, ene, chg)); Move(); Center(); cout << "Energy (" << setw(7) << i + 1 << "): "; cout << fixed << setw(15) << setprecision(5) << Energy() << endl; } } int main(int argc, char** argv) { if (argc < 3) { cerr << "Simple Molecular Dynamics by nox" << endl; cerr << endl; cerr << " Usage: " << argv[0] << " crd_file num_step" << endl; exit(1); } task_scheduler_init init; stringstream ss; int nstep; ss.str(argv[2]); ss >> nstep; SimpleMD md(nstep, argv[1]); md.PrintCrd(); md.Run(); md.PrintCrd(); return 0; }

コンパイル方法

上記のコード中、赤字部分(ピンクを含む)は削除すること。TBBライブラリへのパスを通した後で、Linux(GCCやIntel C++)では、

g++ simple_md_tbb.cpp -O3 -ltbb

icpc simple_md_tbb.cpp -O3 -ltbb

Windows(VC8/9)では、

cl simple_md_tbb.cpp /EHsc /O2 /MD

と、すればよい。

続きを読む...

2009年6月11日木曜日

プログラミングができるということは、人生を楽にできるということ

このブログ記事をはてなブックマークに追加

どのようにプログラマになったのか知るのは面白く、興味深い。そして、プログラミングに対するの様々な考え方や捉え方を知り、それを自分の知識に加えるのが楽しい。そこで自分自身がどのようにしてプログラミングを学んできたのか、そしてプログラミングについてどのように考えているのか述べたいと思う。たいした内容でもないが書いているうちに長文になってしまった。取り敢えず概要だけを知りたい方は、それぞれの段落の最初の文章を読めば何となく分かると思う。

最初のプログラミング

初めて触れたPCは、父親が購入したNECのPC-8801mkIIだった。当時のPCには最初からBASIC(N88-BASIC)が付属しており、これを使って手軽にプログラミングができた。ただ、BASICはインタプリタであり、当時のPCの性能と相まって処理速度が非常に遅く、高速化のためにはアセンブラが必要だった。そこで、父親の書籍を漁りつつ、ニーモニックをハンドアセンブルでマシン語に変換したりもした。

しかしながら、BASICのプログラミングですら当時の自分にとってはとても難しく感じた。その頃読んでいたマイコンBASICマガジンに載っていたプログラムのようなコンパクトでエレガントなコードに比べ、自分のコードはなんて拙くて汚いのだろうと何度も思ったものだ。結局、プログラミングはセンスがある人だけのもので、自分には無理なのだろうかとさえ考えた。あるアイデアがあっても、それを思ったようにコードに落とせないのだ。毎回リファレンスとのにらめっこになる。それでも何とか書き上げたコードは不必要に肥大であり、分かりにくく、あちこちからバグが顔を出していた。

それでも、不細工だろうが何だろうが、プログラムが完成するのはとても嬉しいことだった。自分の力でゼロから何かを作り上げるという行為は本当に楽しかった。後に気が付いたことだが、プログラミングはセンスなんかよりも、この「楽しい」という気持ちの方がよほど重要だったのだ。センスがあっても楽しくなければ長くは続かないだろうし、楽しいと思っているのならば小さな積み重ねが経験となり、知識となっていく。

ただ、その当時はプログラミングばかりをしていたわけではなかった。PCゲームにもはまっていた。オールマシン語が売り文句の一つだった頃だ。そして、その頃のゲームソフトは個人もしくは少数の突出したプログラマによって作成されることが多く、ハイドライドの内藤さんやザナドゥの木屋さんたちのようなカリスマゲームプログラマが存在した。雑誌には特集が組まれ、キャラクタの重ね合わせ処理とか他機種への移植だとかについて語られていたなぁ。このような背景もあって、プログラマには少なからず憧れを抱いていたが、PCを使用できる時間が一日でせいぜい1~2時間と制限されており、その時間の大部分をゲームに充ててしまっていたので、プログラミングにはあまり時間を割いていなかった。

C言語との出会い

大学に入ってからは、自分でPCを持つようになったでの、いつでもPCをいじれるようになり、再びプログラミングを行うようになった。そしてN88-BASICの次に出会ったのがQuickBASICだった。確か雑誌の付録として付いてきたものだ。それはN88-BASICとは一線を画しており、構造化に対応して見た目も大きく異なった言語だった。使ってみるとこれがとても使いやすく、構造化によってプログラムの組み立ても非常に楽になった。これさえあれば、他のプログラミング言語などいらないのではないかと思ったぐらいだ。

その頃、C言語のこともよく耳にしていたのだが、BASICは簡単だし、使いやすいし、すぐに実行できる。そんな素晴らしい言語があるのに何だか面倒くさそうなC言語なんていらないじゃないかと思っていた。特にQuickBASICは素晴らしいと。しかし、すぐにその考えは誤りだと気付かされることになる。

確かこれも雑誌の付録だったと思うのだが、QuickC(Turbo Cだったかな?)に触れる機会を得た。そこまで皆が良いと言うC言語がどれほどのものか少しばかりいじることにした。そして、その簡潔で小気味良い言語設計にすっかり魅了されてしまった。コンパイルしてできあがった実行ファイルはコンパクトで処理が高速だった。これ以降、C言語をメインで使うことになる。

ところで、この時に使った参考書が良くなかったため、C言語の習得が遅れたと思う。特にポインタ関連についてだ。その後、カーニハンとリッチーの書いた「プログラミング言語C」(通称K&R)を読むことで、やっと理解ができるようになった。「K&R」は取っつきは悪いが、実はそれほど分量が多いわけではなく、じっくり読んでみると必要なことを不足なく学べる良書だった。さすがC言語のバイブルと呼ばれるだけはある。これで、C言語を一通り使えるようになったら、ウェブでも書籍でも良いので「C言語FAQ」を読むのが良いC言語の勉強法だと思う。自分は「K&R」と「C言語FAQ」の二冊で間に合ってしまったため他の初学者用参考書についてはそれほど詳しくはないが、一つ言えるのは悪書にだけは手を出さないようにすることだ。最近では評判を聞くに「明解C言語」あたりが良いらしい。まあ、現在ではちょっとネットで調べれば簡単に情報が出てくるだろうからしっかり確認して欲しい。

Perlという新しい扉

さて、この頃にもう一つ、自分のプログラミングに対する考え方に重大な影響を与える出会いがあった。それはPerlだ。Perlであればこれまででは考えられないような手軽さで様々な処理が可能になると聞き、半信半疑ながらも苦労して環境を整えて使ってみた。そして、初めてPerlを使ったときは本当に衝撃的だった。データやファイルを読ませるのも、加工するのもお手の物、なにより環境に密接に結びついており、自分のやりたいことをすぐに表現できることがまさに魔法であった。これはラクダ本の日本語訳が出る前のことだ。そして、Perlの虜になった自分はC言語とともに長く使い続けることになる。

ところで、新しい言語に出会うたびに自分の視野の狭さを痛感する。それまで知らなかった世界がいきなり目の前で開かれるような感じだ。自分は今までなんて狭い世界に住んでいたのだろうと。無理やり喩えると、美味しい料理と言えばステーキだったが、寿司も旨いことを知るようなもの…かな?

オブジェクト指向の世界: C++、そしてPythonへ

C言語とPerlで自分のスタイルはほぼ完成したかに思えたが、技術革新に終わりはなかった。その後の技術でもっとも大きく影響を受けたのはオブジェクト指向プログラミングだ。これは、部品であるオブジェクトでプログラムが構成され、クラスによる継承、カプセル化、多態性などの機能を有する。部品を構成して作るというイメージがプログラミングによりモノを作るという考えにピッタリ一致しているのが個人的にはツボだった。

ここでC++が登場する。これは、C言語を拡張した言語として1980年代に生まれた。C言語の仕様はそのままにオブジェクト指向プログラミングのための機能を追加した言語となっている。自分が本格的に使い始めたのはVisual C++ 1.0からだったが、最初の頃は主にベターCとして使っていた。というのもオブジェクト指向プログラミングについてちゃんとは理解しておらず、さらにC++自体がオブジェクト指向プログラミングできる言語であったため、オブジェクト指向プログラミングだけできる言語とは違い、オブジェクト指向の知識のないままコーディングをすると簡単に構造化プログラミングとごちゃ混ぜになってしまうからだ。案の定、見よう見まねで作ったオブジェクト指向まがいのコードはそれはひどいものだった。

それでも、オブジェクト指向のコードを読んだり書いているうちにそれなりに解るようになり、特に「デザインパターン」によって理解が深まった。デザインパターンで有名なのはGoF(Gang of Four)の「オブジェクト指向における再利用のためのデザインパターン」だろう。これには本当に感動した。オブジェクト指向の持つ本当の力をデザインパターンで知ったと云っても過言ではない。こうして、自分の中にオブジェクト指向の考えが浸透していった。

一方、これまで使ってきたPerlであるが、不満もあった。読みにくいのだ。Perlハッカーからすればそんなことはないのだろうが、それほど使いこなしているわけではない自分では今まで書いてきたコードでさえ読みにくかった。初期の頃、しばらくたった頃、最近…と、時期によってスタイルが変わっていた。もちろん書き方の自由度が高いのはPerlらしさでもあるのだが、自分にとっては少々鬱陶しく感じてしまっていたのだ。あと、記号を多用するのがあまり好きではないこともあった。

オブジェクト指向プログラミングを知ってしまってからは、Perlに替わるオブジェクト指向の使える良いスクリプト言語がないだろうかと考えるようになった。そこで目を付けたのがPythonだ(Rubyという選択肢もあったがそうはしなかった。理由については「ビューティフルコードとmalloc/free論争」に書いてあるので興味のある方は読んで欲しい)。もともとPythonはバージョン1の頃に少し使ったことがあったのだが、そのときは慣れないインデントの仕様もあってそのまま使わなくなってしまった。本格的に使い始めたのはバージョン2.3以降である。初めのうちは違和感を覚えたインデントも慣れていくにしたがって、もうこれしかないとまで思うようになった。というのも、自分の求めていた読みやすさに大きく貢献するスタイルだったからだ。誰が書いても同じような見た目になるのは素晴らしい。誰が書いてもというのは過去の自分でもということだ。自分で書いた使い捨てレベルのコードでも何を書きたかったのか明瞭だった。そして、クラスを使ったオブジェクト指向プログラミングも完備していたので、少し大きめのコードを書くときにも困ることはなかった。すぐに、ほとんどの雑多な仕事をPythonに任せるようになった。Google App Engineでも使えるしね。

こうして、CからC++へ、PerlからPythonへと変わっていったのだ。

ところで、C++を覚えるためには、ある程度の覚悟を持って学ぶ必要があるかもしれない。C言語を知っているのであれば、初めはベターCのような使い方でも良いと思う。その後、クラスに慣れていき、オブジェクト指向について学ぶのが良いのではないだろうか。オブジェクト指向についてよく知っているのであればそれほど難しくはないだろう。ただし、C言語に関する知識は必要となる。また、STLの理解のためには「STL標準講座」が役立つ。ある程度C++について理解したのなら、まず「Effective C++」を読もう。そして前述のGoFの「オブジェクト指向における再利用のためのデザインパターン」、さらに理解を深めるために「Modern C++ Design」を読むといいだろう。

Pythonについては勉強らしいことはほとんどしなかった。せいぜい公式サイトのドキュメントにあるチュートリアルぐらいである。数日もあれば十分ではないだろうか。あとはリファレンスさえあれば困ることはない。とは言っても、深く理解するにはそれなりに使い、学ぶ必要はあるけど。書籍では「Pythonクックブック」があると便利だと思う。

プログラミングは目的ですか? 手段ですか?

さて、こうして自分は日常的にプログラムを書くようになった。では何のためにプログラムを書いているのか? プログラミング自体が楽しいから? プログラムを作る必要があるから? それとも他の理由から? そもそも、プログラミングは目的なのか? 手段なのか?

趣味であればプログラミングは目的でも手段でもあると言える。それらすべてを含めて楽しいわけだ。ただ、仕事やもっと大きな視野からみた場合、プログラミングは手段でありたいと思う。先に大学に入ったと述べたが、自分の専門はコンピュータ関連ではなく薬学だ。そして、薬剤師でもある。薬学を選んだのは自分がもともと化学に興味があり、人々に直接役に立つ化学が薬学だと考えたからだ。そして、将来的に薬学をコンピュータによってさらに発展させられるのではないかと感じていたからだ。現在は薬学系の研究の面からプログラミングが役に立っている。仮に研究ではなく調剤をしていたとしても、プログラミングにより多くの恩恵を受けることはできただろう。

今の世の中、多くのアプリケーションソフトウェアがあり、わざわざプログラミングする必要なんてない、プロではないのなら意味がないとも言われたりするが、本当にそうだろうか。自分自身をみると仕事でもプライベートでもプログラミングは非常に役に立っているし、道具(手段)としてとても重宝している。仕事にしてもプログラミングあっての仕事ではなく、プログラミングによって便利で楽に素早く進めることができるというだけだ。喩えると、重い荷物を届けるのに歩いて(=コンピュータを使わない)運ぶか、電車やバス(=アプリケーションソフトウェア)で運ぶか、自家用車(=プログラミング)で運ぶかという違いに似ている。自家用車であれば自由にどんなところにも運べるわけだ。今、社会はコンピュータと密接に繋がっている。そして、それらのコンピュータと自由に意思疎通するためにプログラミングがある。

自分にとってプログラミングができるということは、人生を楽にできるということなんだと思う。

続きを読む...