TELEMAC parallel

TELEMAC のインストール時に、mpich2 が含まれていました。

「OpenMP ではなくて MPI か」で素通りしました。というのも、MPI の概要は理解していましたが、具体的な利用手順、コマンドを一切知りませんでした。

せっかくなので telemac2D を parallel で動かそうと思ったのですが、そのコマンドがわかりません。マニュアルを見ても具体例が書かれていません。で、Forum を検索。発見。

まずは初回にユーザーと pass を登録。

mpiexec -register

次に実行。

CONFIGNAME '-c wing64mpi' を加えただけです。CONFIGNAME って、config ファイルに書かれているコンパイラー識別コードのようなモノのことだったのですね。 

telemac2d.py -c wing64mpi ***.cas

cas ファイル内で4コア指定

PARALLEL PROCESSORS = 4

これで parallel 版も動きました。single で14.5時間かかった自由表面流れのモデルが、4コアで9.5時間でした。早くなっています（といっても、２次元計算にしては遅いと感じています）。

Forum を見る限り、領域分割を行っているようです。今、詳細を追いかけようとは思っていませんが、いずれ理解したいですね。

R apply の問題

Rでは、for 文の代わりに apply 群を使われる方が多いようです。

確かに、apply 群を使って書く方が、簡素で読み易いと思います。また、for 文より速くなる場合があるようで、web上ではしきりにお勧めされていました。（apply の中で for が使用されていましたので、劇的に、ということはないと思いますが）。

念のため、apply や lapply を使用して書き直し。
もともと、dependency を算出する関数が、ベクトルを扱えなかったこと、並列化をもくろんでいたことからfor を使っていました。今回は、関数の引数にリストから値を代入する関数を作り、それを apply や lapply で呼んでみました。

で、計算!

が、ダメ。

夜中に計算をかけて、朝に「終わったかな？」と確認しようとすると、動きが異常に遅い。SSDにスワップファイルを作っているようです。リソースモニターを見ると、物理メモリ24GBを使い切っていました。for や foreachであれば、ループのたびに変数を入れ替えながら計算が進むため、それほど大きなメモリは必要としません。が、apply 群では大きなマトリックスをそのまま保持して計算が進んでいくため、いくつかの計算ステップが進めば簡単にメモリ不足に陥るのでしょう。
計算が進むたびに途中経過をディスクに書き出せばメモリを解放できますが、遅くなります。ま、扱うデータの容量を考えながら、方法を選択しないといけないのでしょうね。ビッグデータを扱う場合のパッケージもあるようですので、必要になった段階で検討してみましょう。

いずれにしても、今回は並列化のほうが速いという結果になりました。

コア数　　　計算時間　メモリ
for（Single）　11.5h　　300MB
foreach（4）　　3.7h　　1.5GB
foreach（6）　　3.8h　　1.8GB
apply（Single）6hで中止 20GB以上

これで進めてみましょう。

R foreach で並列化

R で並列化・高速化を扱う場合、以下の2通りがあるようです。

1. 関数をC/C++ に移植。さらにCUDAを利用
2. foreach による並列化

頻繁に使うツールであれば前者でしょうが、今回はそこまでの頻度・意欲がありません。で、2を選択。

最初はビルド時の自動並列化のようなオプションがあるのか？と思い探してみましたが、見当たりません。代わりに見つけたのが 2 の foreach 。
R では、for を foreach に変更するだけで並列化してくれるようです。容易なためか、メジャーな手法のようですね。

先日のコードに対し foreach を使って並列化。
が、今度は排他制御に関するコマンドが見つかりません。複数のスレッドから同一変数への書き込みを避け、ファイルに追記するよう変更したのですが、タイミングが重なった場合にデータが飛んでしまいます。このあたりを制御するコマンドが見当たりません。また、foreach 内の break も効かないようです。

そのままではうまく動いてくれなかったので、部分的にコードを書き直し。

っで、計算！

今度はうまくいきまた。

1 コアで 11.5 時間だったのが、4 コアで 3.7 時間まで短縮できました。
6 コアでも計算してみましたが、頭打ち。メモリへのアクセスが集中しますので、遅くなったのでしょうか。
ま、それほど早い計算時間でもないのですが、このあたりが落としどころでしょうね。これ以上は 1 の方法、R から C へ移植しするしかないでしょう（実務なら、迷わずそちらを選びますが）。

扱いたいデータがもっと大きなサイズですので、3.7時間という結果は大きな制約です。データの前処理を考えないとだめでしょうね。なんだかベクトル化の方が早いかも？と思えてきました。

で、もう一度チャレンジ。

今度はクリアーできましたが、新たに問題発生。


続きは後日。

スタックサイズの変更

並列化済みの Dtransu を動かす場合、スレッド数やスタックサイズは、batファイル内で指定していました。

後輩から、それをソース内で固定できないか？といった問い合わせがありました。
ソース内で固定したことはないのですが、当然可能だろうと触ってみたところ、うまくできません。力不足です。

調べてみると、スレッド数の指定はすぐに引っかかりました。Win7 + Fortran では以下の通り。

call OMP_SET_NUM_THREADS(スレッド数)


しかし、スタックサイズの指定が分かりません。オプションの中を見てもわからなかったので、その日は解決できず。


後日、調べてみましたが、結局明示する方法はわからず。最終的に、Win7側で環境変数を固定することにしました。

OMP_STACKSIZE　を追加すればOKです。スタックサイズはそれほど変更することがないので、固定しておいても問題ないですね。画像では、OMP_NUM_THREADSも追加していますが、ソース内で指定する場合は不要です。

しかし、これらを解決するのに、3夜かかりました。
プロがそばに欲しいですね。

アプリケーションを正しく起動できませんでした(0xc000007b)。

以前、Dtransu の OpenMP 64bit版を他の PC で動かした時のエラーについて、書き残していました。
http://phreeqc.blogspot.jp/2011/11/libiomp5mddll-openmp.html
 
今日、その流れで作業していた後輩から一言、「動きません。エラーでます」。
 
 うーん、と思いつつ、エラー内容を聞いていると、やはりdll。libiomp5md.dll は問題ないので、それ以外の dll に関するエラーです。1つずつ調べようかと思いましたが、ちょっと込み入っていましたので、以前教えていただいた「Redistributable Libraries for 32-bit/64-bit (x64) msi files」を DL して渡しました。

これをインストールすると、問題は全て修正され動くようになったようです。
最初から、こうすれば良かったですね。

libiomp5md.dll と OpenMP

急ぎの計算があったので、後輩の Workstation を借りています。

Del l の precision T3500です。
残念ながら4コアでした。

早速、OpenMP で並列化済みの Dtransu とデータファイルをHDにコピーし、起動！

と思いきや、いきなりエラーが出ました。

「コンピューターにlibiomp5md.dllがないため、プログラムを開始できません。・・・」

そういえば、コンパイル時にライブラリーにリンク張る設定にしてました。（っていうか、並列化コード、使ってないんですね。もったいない。）

libiomp5md.dll を自分のpcからコピーして再度実行すると、以下のエラーが出ます。

「アプリケーションを正しく起動できませんでした(0xc000007b)。・・・」

うーん、これは前にも見おぼえがあります。でも、忘れました。

検索してみると、X86の dll を64bit OS で使用した場合に出るエラーでした。

64bit 用の dll に入れ直すと、正常作動しました。

うーん、OpenMPで並列化したコードを配布する場合、このライブラリーは一緒に配布して良いのでしょうか？サポートに問い合わせると、以下の回答が得られました。

「可能リスト(fredist.txt)に入ってますし、Redistributable Libraries に入っているので問題ないですよ。」

そうですか。そんなリストがあったんですね。
確かに、書かれていました。どこまで見ていないんでしょうか。



余談ですが、SSRの更新をしたにもかかわらず、登録をしていなかったため、Update 7もDLできない状態でした。あわてて登録、VFを12に、shellを2010に入れ替えました。はあ、プロが側に欲しい。

さらに余談ですが、12月の地下水学会の講習会で2次元版DtransuのOpenMP並列化コードが配布されるようです。他支店から参加するので、参考に見せてもらいましょう。


------------------------------------------
2012.8.24追記
上記だけではダメな場合もあります。
基本は Redistributable Libraries をインストールしてしまいましょう。
http://phreeqc.blogspot.jp/2012/08/0xc000007b.html

今年の目標

今年も半分終わりました。

短期目標として、今年は以下の点をこなそうと思っていました。進捗は50%程度と言ったところでしょうか？

1. Dtransuの高速化　OpenMP, GPGPU
2. 飽和・不飽和輸送と反応の連成
3. 変形解析

Dtransuは半分クリアーです。あとはGPGPUだけですね。
ただし、Dtransuのコード読みが終わっていないので、CUDAコードと併せて読み込む必要がありそうです。目標にはしていませんが、使用している技術者として当然のことなので、やっておく必要があります。でも、難しいんですよね。

輸送と反応は1次元のみHP1でクリアーです。3次元はMATLAB+COMSOL+IPhreeqc と連携コードで可能ですが、ソフトがありません。COMSOLのセミナーには出席する予定ですが、ペイできなければ買ってくれないでしょう。
しかし、お金を出せばできることは分かりました。この辺りでクリアーにしましょう。

変形解析は、まだまだですね。
分かったと思ったら、次に分からないことがドンと出てくる段階です。山の裾野のほうでしょう。早く連続体から不連続体、粒子法へ行きたいのですが、先は長そうです。学ぶべきことは、まだまだたくさんあります。

残り半年、焦らずに時間をかけてやりましょう。

cloudPEST

USGSがユニークなモジュールを公開しました。

cloudPESTです。
http://pubs.usgs.gov/of/2011/1062/
Amazon Elastic Compute Cloud (Amazon EC2)を使って、PESTを動かそうという試みのようです。
私はEC2を使用したことはありませんが、並列計算でコアを追加したい場合には使えるかもしれません。

今日は樫山先生の「並列計算法入門」を読んでいました。スペックについては少し古い感はありましたが、理論は今も変わっていません。というより、実務では追いついてないレベルでしょう。陽解法、陰解法での領域分割の仕方、差分法やFEMでの並列化の考え方など、今まで「もやっ」としていた頭が整理できました。
並列化を目指すには、MKLやHEC-MWなど、既に出ているモジュールを利用していくのが実用的なのでしょう。（まだGPUの使用を目指してますが。）

クラウドや並列化、環境は整っていますが、なかなか先に進みません。

最先端の数値解析

地盤工学会誌3月号が届きました。

「最先端の数値解析」といった興味を引く内容で、すぐに開けてみました。
まず目を引いたのが、冒頭の若井准教授の「無理して”結果”をあわせないこと」。思わずドキッとして呼んでいくと、なんとなく感じていたことが明確になりました。最後の「3.完全に合わなくて良い」が個人的に好きですね。甘く感じる章題に反し、「実測に合う予測結果が得られる組み合わせが”全く”ないのかどうかが、個々の数値解析法の良否を判断する上で重要な材料となりうる。」と非常に厳しい正論を掲げられています。つい、契約金額と時間のバランスに逃げてしまいがちですが・・・反省ですね。

あとは同じ群馬大の「マルチコアに対応する浸透流解析・・・・」ですね。今日、偶然ですがソフト会社さん2社より電話をもらい、並列化やGPUコーディングの話をしたところでした。業界ではプチブームなのでしょうか。過去の資産を生かすには越えないといけない壁でしょうから、どの会社も取り組むべき課題なのでしょう。64bit対応と共に。
ちなみに、このMKL、手元にあるんですが並列化の内容を全く見ていませんでした。今後何かの役に立つかもしれません。明日、確認してみましょう。

領域分割による並列化

ダイヤコンサルタントさんのHPを見ますと、並列化したDtransuにより600万要素の浸透流解析の実績があるとのこと。ただ、移流分散については書かれていませんね。
http://www.diaconsult.co.jp/leaflets/pdf/dia_k_3Dshintou2.pdf

こちらに記載のあるように、HEC-MWの利用で実現しているのでしょう。
http://www.ssken.gr.jp/MAINSITE/download/newsletter/2009/20090903-sci-1/lecture-5/ppt.pdf
やはり、領域分割ですね。

前にも書きましたが、領域分割による並列化は樫山先生がお得意の分野です。収束計算ではここがメインにならざるを得ないのでしょう。

HEC-MWはフリーですので、環境としては整っているわけですが、そちらに移植するかどうかというと迷います。どちらかというと、Fortranのまま、領域分割のアルゴを追加し、GPU対応コンパイラーでビルドするという方針のほうが、将来性があるように思えます。同時にOpenMP対応版もできそうですからね。ま、どちらにしても、プログラマーが必要でしょう。

海外だと、すでに並列化済みのものが公開、あるいは販売されているかもしれません。飽和、差分ではVisual MODFLOW (MT3DMS)が並列化済みのようですので、不飽和帯との分割計算も手かもしれません。
http://www.swstechnology.com/groundwater-technical-support/software-faqs-prod?ProdID=29

探してみましょうか。

GPGPU・CUDA

GPUで並列計算、すごそうですよね。

計算中遊んでいるGPUを活用できるなら、劇的に早くなるでしょう。CPUを使った並列計算ではコア数や今後の開発に限界がありますが、GPUだと並列処理能力が違いますし、現在でも伸びている分野です。CPU交換するよりも簡単ですし、期待しちゃいますよね。GPU対応のfortranコンパイラーも販売されています。
http://www.softek.co.jp/SPG/Pgi/Accel/
OpenMPと同じよな感覚でコマンドをはさんでいくだけのようです。それだけでGPUが使えるなら便利ですよね。

でも、結局はどこを並列化するか適切に判断しないと、パフォーマンス向上は期待できないんですよね。特に収束計算では難しい。

うーん。樫山先生がされているように、計算領域の分解で並列化するしかないのかな？やり方詳しくわかりません。土壌物理が終わったら、ちょっと調べてみましょうかね。

Dtransuの並列化！

コンパイラーの最新版がきたので、とりあえず、ビルド！
今組んでいる23万点モデルの初期定常で計算時間をチェックしてみました。

まず、/O2を/O3で試しましたが、速度が上がらず。次にSSE4.2を入れると微妙に早くなりますが、殆ど変化しない状況でした。これには期待していたので残念です。自動並列化に関しては逆に速度が落ちてしまいました。結局、これらの設定を変えただけでは速度の上昇が認められませんでした。

繰り返し計算では、収束したら次のステップへという設計のため、loopの回数が決まっていないことも並列化しづらい点であるように思われます。最終はOpenMPですが、ここまでくると時間が掛かるため、ソフト屋さんにお金を払ってでもやってもらいたいですね。明日以降も、色々触ってみようと思います。