ROOT で OpenMP を使う

最近は dual corequad core の時代になったのに、自分でその性能を引き出すのは困難です。CPU の core 数ばっかり増えても、解析の速度は大して速くなっていない。というわけで、OpenMP を使って ROOT でも並列処理をできないか、というお話。

OpenMP とは

OpenMP というのは、難しいこと考えてなくても、あなたの code を裏で multi thread 化し、multi core の CPU の性能を引き出してくれる便利な道具です。細かいことは検索して調べて下さい。

OpenMP と Minuit2

ROOT の解析で時間を食うのが、多くの場合 fitting でしょう。あと思いつくのは、巨大な tree や ntuple から数値を読み出すとき。前者は CPU 速度に比例しますが、後者は多分 disk I/O の方が効くでしょう。後者は SSD にするなりが必要ですが、前者は multi core を活かすことで高速化できます。

さて、ROOT だと Miniut という library が標準的に使われます。こいつの新しくなった Minuit2 というのがあって、特に目新しい部分は、OpenMP に対応しているというところです。Minuit と Minuit2 の他の差を理解していませんが、OpenMP 関連のところに関しては、多次元空間で行われる fitting の minimization を、並列化して高速化してやろう、ということだと思います。

さて、OpenMP に対応した Minuit2 を ROOT で使うためには、次の条件が必要です。

  • GCC 4.2 以降を使うこと。GCC 4.2 以降であれば、OpenMP に標準で対応しているため。
  • ROOT を build するときに、Minuit2 と OpenMP を有効にすること*1
  • Fitting の engine を手動で Minuit2 に切り替えること。

まず、ROOT を configure するときに、--enable-minuit2 というのを付けましょう。これで Minuit2 も build されます。ROOT の compile については、ここを参照して下さい。

$ ./configure --enable-minuit2

次に、make を始める前に、以下の環境変数を設定しておきましょう*2

$ export USE_PARALLEL_MINUIT2=1
$ export USE_OPENMP=1

最後に、make します。上記環境変数を設定しておくと、g++ の引数に "-fopenmp" がついたりしてくれます*3。既に build 済みの ROOT がある場合は、$ROOTSYS/math/mathcore/src/FitUtilParallel.?、$ROOTSYS/math/minuit2/src/Numerical2PGradientCalculator.?、$ROOTSYS/lib/libMathCore.*、$ROOTSYS/lib/libMinuit2.* あたりを消去してから make し直しましょう。

さて、この OpenMP を使って build された Minuit2 を使ってみましょう。

#include "TH1D.h"
#include <iostream>

Double_t fit()
{
  TH1D* h = new TH1D("h", "h", 100000, -5, 5);
  h->FillRandom("gaus", 10000000);
  Int_t start = clock();
  h->Fit("gaus", "l");
  Int_t stop = clock();

  delete h;

  Double_t t = (stop - start)/1e6;
  std::cout << t << " (sec)" << std::endl;
  return t;
}

この script を Minuit と Minuit2 で実行してみます。下記の例は、Minuit2 を使った場合です。TVirtualFitter::SetDefaultFitter から Minuit2 を指定します。

$ root
root [0] .L fit.C+O
root [1] TVirtualFitter::SetDefaultFitter("Minuit2")
root [2] TH1D* h = new TH1D("", "", 10, 0, 0)
root [3] for(int i = 0; i < 100; i++) h->Fill(fit());           
root [4] h->GetMean()

結果は、Minuit2 が平均 0.740 秒、Minuit が平均 0.656 秒となりました。あれ、Minuit のほうが速い。Minuit2 でやると、確かに CPU 使用率は 120% くらいに上昇します。ただ、実質的には遅くなってしまいました。もうちょっと複雑な fitting や error の評価をやると効いてくるかもしれません。要検証*4

ACLiC から OpenMP を使う

さて、残念ながら Minuit2 は期待した結果になんなかったですが、次は ACLiC で OpenMP を使う方法です。林檎生活100: OpenMP 入門 #3 を参考にして、ROOT の ACLiC で π を計算してみます。

次の code を、pi.C という名前で保存します。

#ifdef private
#undef private
#endif

#include <omp.h>
#include <time.h>
#include <iostream>

void pi()
{
  const int kN = 100000000;
  long double top = 0.L;
  long double bottom = 0.L;
  long double kStep = 1.L/kN;
  long double sum = 0.L;
  
#pragma omp parallel
  std::cout << "Hello, OpenMP!\n";
  
#pragma omp parallel private(top, bottom)
  {
#pragma omp for reduction(+:sum)
    for(int i = 0; i < kN; i++) {
      top = 4.L/(1.L + (kStep*i)*(kStep*i));
      bottom = 4.L/(1.L + (kStep*(i+1))*(kStep*(i + 1)));
      sum += (top + bottom)*kStep/2.L;
    } // i
  }
  std::cout.precision(16);
  std::cout << "pi = " << sum << "\n";
}

ROOT を起動し、次の内容を打ちましょう。

$ root
root [0] TString cmd(gSystem->GetMakeSharedLib())
root [1] cmd.ReplaceAll("g++", "g++ -fopenmp")
root [2] gSystem->SetMakeSharedLib(cmd)
root [3] .L pi.C+
root [4] pi()

0 行目は、ACLiC 経由で script を compile するときの、default の引数を取得しています。GCCOpenMP を使うときには -fopenmp という引数が必要なため、1 行目で "g++" を "g++ -fopenmp" に置き換えています。そして、置換した結果に再設定するのが 2 行目です。3 行目で pi.C を compile しつつ load し、4 行目で関数 pi() を呼び出しています。

そうすると、このように π の計算結果が現れるはずです。

Hello, OpenMP!
Hello, OpenMP!
pi = 3.141592653589793

"Hello" が 2 度現れているのは、2 つの CPU core で 2 つの thread が走っているからです。1 CPU で計算させるよりも、およそ 2 倍程度、計算が速くなっているはずです。pi.C の中の、"#pragma" の行を全て comment out して試してみて下さい。

ちなみに、pi.C の先頭の 3 行は、ACLiC で実行するときに必要になります。private という予約語がどっかで使われてしまうため、それを #undef しています。ROOT を使わず、普通に C++ だけでやる場合には不要な行です。

*1:Binary distribution がどうなってるのかは知りません

*2:bsh や zsh の場合です。tcshcsh のときは、setenv して下さい。

*3:OpenMP を使うときに GCC で必要な引数です。

*4:こないだ pol4 でやったときは速くなった気がしたんだけど。