最近はマルチコア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倍
さすがに複数コアを利用するとある程度はサチってしまうが、それでもこれだけ簡単に利用できることを考えると十分な性能が発揮されていると思う。
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
と、すればよい。
これまでは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
最後に実際の計算時間を測定したので、それを示す。チャートではコアに対する実行速度向上の割合を示している。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倍
さすがに複数コアを利用するとある程度はサチってしまうが、それでもこれだけ簡単に利用できることを考えると十分な性能が発揮されていると思う。
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
と、すればよい。
コメント