計算工学ナビ

ものづくりにHPCを活用するための ツールとケーススタディー

サイト内検索

Windows10のubuntu互換環境でオープンソース流体ソフトウェアを動かしてみた

東京大学生産技術研究所 革新的シミュレーション研究センター

鵜沢 憲 (uzawa@iis.u-tokyo.ac.jp)

14274439_10208515693423093_1321850254_o

背景

2016年8月3日に提供された64ビット版Windows10のバージョン1607(ビルド番号14393, Anniversary Update)から,Windows上でubuntu互換の実行環境(Windows Subsystem for Linux, 以降WSL)[1],及びシェル(Bash on Ubuntu on Windows, 以降BUW)が使用可能となった.

これまでWindows上で流体ソフトウェアを動かす場合,Windowsのバイナリで提供してもらって実行する,もしくはVMware[2], VirtualBox[3]等の仮想環境やCygwin[4]等のエミュレータでUNIXのバイナリを実行するようなケースがほとんどであろう.

Microsoft社のWSL開発グループによれば,それらの仮想環境と比較して,今回リリースされたWSLは必要とするリソース(CPU, メモリ, ストレージ)が少なくて済むと述べている[5].また,これまでの仮想環境やエミュレータで不便を感じることもあったWindowsとLinuxとの間でのファイルやディレクトリ,コマンドの互換性の問題も回避されている.これらの特徴は,Windowsで流体解析を行う際のハードルが下がる可能性がある.

そこで,試しにオープンソースの流体ソフトウェアをインストールし,流体解析の一連の流れを通じてWSLの使い勝手を調べるともに,簡単なベンチマーク問題を対象に計算性能を評価し,WSLにおける流体解析の可能性を調査することとした.

なお,流体ソフトは文科省フラッグシップ2020プロジェクト重点課題8「近未来型ものづくりを先導する革新的設計・製造プロセスの開発」サブ課題A「設計を革新する多目的設計探査・高速計算技術の研究開発」で開発が進められているFrontFlow/violet-Cartesian(FFV-C)[6]を用いた.FFV-Cは任意のアーキテクチャ対応の性能評価ライブラリPMlib[7]をデフォルトで実装しており,計算終了後に計算時間の測定や計算負荷のホットスポット同定等を検証することができる.

今回,PMlibを利用してWSLでの計算時間や計算性能について調査したところ,いくつか興味深い結果が得られたので,WSLの使用感とともに,まずは第一報としてここに報告する次第である.

WSL及びBUWのインストール及び動作確認

WindowsへのWSL及びBUWの導入手順や不具合,使用感などはこれまでに多く報告されている[8-11].

今回筆者も上記の手順をなぞることにしたが,日頃の行いが悪いのだろう,第一歩目のバージョン1607へのアップデートに失敗してしまった.プリインストールされていたWindows10 HomeにAnniversary Updateをかける度に,途中でフリーズしてしまうのである.頑張って心折れずに調べたところ,システムをSSDに入れているPCにアップデートをかけるとフリーズする場合がある[12] とのことで,今回の原因もどうやらそのようであった.

対処方法として,今回はプリインストールされていたWindows10を完全に削除し,イメージ形式で取得したバージョン1607[13]をUSBからクリーンインストールすることとした.クリーンインストール後は,上記の記事通りになぞることで,無事WSL及びBUWを入手することができた.

インストールが上手くいくと,スタートメニューに「Bash on Ubuntu on Windows」が現れるようになる(図1).インストール時間を正確に測定したわけではないが,感覚的には,従来の仮想環境やエミュレータ環境を構築する場合と比較してはるかに短時間かつ容易に構築することができた.一般的なWindowsアップデートに毛が生えた程度の印象である.

図1Bash on Ubuntu on Windows on スタートメニュー

図1Bash on Ubuntu on Windows on スタートメニュー

クリックすると,普段使っているubuntuと同様のbashが立ち上がった.ubuntuとbashのバージョンは,それぞれ14.04と4.3.11である.(図2)

図2 夢にまで見たWindowsでのbash

図2 夢にまで見たWindowsでのbash

BUWを一通り触り,ネイティブbashと同様の操作感であることを確認.

FFV-Cのインストール及び動作確認

引き続き,流体ソフトウェアのFFV-Cをインストールする.なお,WSLには必要最低限の機能しか入っていないので,FFV-Cの計算やプリ・ポスト処理に必要なソフトウェアやライブラリを予めインストールしておく.

GNUコンパイラ(gcc, g++, gfortran)及びgnuplotのインストール

ネイティブubuntu下と同様に,それぞれsudo apt-get install * とすれば良い.

$ sudo apt-get install build-essential
$ sudo apt-get install gfortran
$ sudo apt-get install gnuplot-x11

(もちろん,まとめてインストールしても良い.)

なお,sudo実行時に「***の名前解決ができません」と出るときは,WindowsのIPアドレスをetc/hostsに記入すれば,毎回怒られずに済む.(図3)

図3 /etc/hosts ファイル

図3 /etc/hosts ファイル

OpenMPIのインストール

FFV-Cは,OpenMPIによるプロセス並列とOpenMPによるスレッド並列を組み合わせたハイブリッド並列でコーディングされている.したがって,OpenMPIのソースコードを公式サイト[14]等から入手し,所望のディレクトリ下で展開,コンパイル・実行しておく.

この際configure用に以下のスクリプトを用意しておくと便利である[15].

 

$ vi config.sh
#!/bin/bash
export CC=gcc
export CFLAGS=-O3
export CXX=g++
export CXXFLAGS=-O3
export F77=gfortran
export FFLAGS=-O3
export FC= gfortran
export FCFLAGS=-O3
./configure --prefix=$1

問題なくconfigureできたら,コンパイルとインストールをする.

 

$ make
$ make install

MobaXterm のインストール

プリ・ポスト処理のために,Windows上の代表的なXサーバであるMobaXterm[16]をインストールする.MobaXtermをインストールしておくと,入力ファイルを外部エディタで修正したり,計算結果を可視化したりできるようになるので便利である.

インストールはほぼ一直線なので割愛する.無事にインストールが完了すると,MobaXtermを立ち上げることができる.(図4)

図4 MobaXterm の初期画面

図4 MobaXterm の初期画面

MobaXterm を確認すると,WindowsのIPアドレスが192.168.103.2:0.0と表示されているので,BUW側で

$ export DISPLAY=192.168.103.2:0.0

と設定する.これでXサーバが利用可能になる.

V-Isio のインストール

ポスト処理のために,軽量のオープンソース可視化ソフトウェアV-Isio[17]もインストールしておく.手順は以下の通り.

 

alien パッケージのインストール

$ sudo apt-get install alien

alien で rpm 形式のパッケージを deb 形式に変換

$ sudo alien Visio-2.4-6.el6.x86_64.rpm

成功するとvisio_2.4-7_amd64.debができる.

 

deb パッケージのインストール

$ sudo dpkg -i visio_2.4-7_amd64.deb

パスを通す

$ export PATH=/usr/local/Vtools/bin:$PATH

事前に足りないライブラリを入れておく

$ sudo apt-get install libgtk2.0-0
$ sudo apt-get install libpangox-1.0-0
$ sudo apt-get install libjpeg62

これでも

$ Visio: error while loading shared libraries: libtiff.so.3:
 cannot open shared object file: No such file or directory

と怒られるので,libtiff.so.* を適当に探し,シンボリックリンクを作成する.

$ sudo find . -name libtiff.so.*
$ ./usr/lib/x86_64-linux-gnu/libtiff.so.5
$ sudo ln -s /usr/lib/x86_64-linux-gnu/libtiff.so.5 \
/usr/lib/x86_64-linux-gnu/libtiff.so.3
$ Visio

アクセス許可を求めるポップアップが出るが,ここで「はい」とすることで,無事にV-Isioが立ち上がる.(図5)

図5 V-Isio の初期画面

図5 V-Isio の初期画面

 

今回は割愛するが,V-Isioの他にも,ParaViewやVisIt等の可視化ソフトウェアもインストールできるので,必要に応じて試して頂きたいと思う.

FFV-Cのインストール

GNU コンパイラ(もしくは Intel コンパイラ)と OpenMPI があれば FFV-C のインストールが可能になる.

$ http://avr-aics-riken.github.io/ffvc_package/

で最新版FFV-C(バージョン2.4.3)をtar.tz形式もしくはzip形式でダウンロードしたのち,所望のディレクトリに展開し,インストール.

なお,FFV-C には GNU 環境用と Intel 環境用の両方に自動インストールスクリプトが付属しているので,これを利用すると良い

$ tar xvfz avr-aics-riken-ffvc_package-x.x.x.tar.gz
$ cd avr-aics-riken-ffvc_package-x.x.x
$ ./install_intel.sh

お茶を一杯飲んでいる間にインストールが終わる.

試しに組み込み例題の2次元キャビティ問題を計算する.プリ処理としてviエディタで入力ファイルに問題ないことを確認した後,実行.

無事にFFV-Cが動作することが確認できた.(図6)

図6 FFV-C on Bash on Ubuntu on Windows

図6 FFV-C on Bash on Ubuntu on Windows

計算ログをgnuplotで確認.(図7)

速度の発散値や相対残差ノルム値の時間発展を確認し,収束性に問題がないかどうかを確認する.

図7 gnuplot を用いた計算ログの表示

図7 gnuplot を用いた計算ログの表示

速度の絶対値をV-Isioで可視化してみる.(図8)

物理的に奇妙な挙動を示していないことを確認する.

図8 2次元キャビティ問題における速度の絶対値のプロット図.ubuntuにインストールしたV-Isioで可視化

図8 2次元キャビティ問題における速度の絶対値のプロット図.ubuntuにインストールしたV-Isioで可視化.

ここまでで,プリからポストまでの一連の流体解析が問題なく実行できることが確認できた.

 

ところで,今回のWSLで筆者が個人的に大変便利だと感じていることは,ubuntuのファイルやディレクトリをWindows操作できる点である.このことにより,一連の流体解析の操作性が従来の仮想環境やエミュレータと比較して格段に向上している.本稿はこのために書いたといっても過言ではない.

WSLは,Windowsのファイルシステムの

C:¥Users¥「username」¥AppData¥Local¥lxss

に格納されている.(図9)

ただし,これを見るためには,「エクスプローラー」→「オプション」→「フォルダーと検索のオプションの変更」→「フォルダーオプション」→「表示」で,「保護されたオペレーティングシステムファイルを表示しない(推奨)」のチェックを外しておく必要がある.(図10)

図9 Windowsファイルシステム内のWSL格納箇所

図9 Windowsファイルシステム内のWSL格納箇所

図10 フォルダーオプション内のチェックを外す

図10 フォルダーオプション内のチェックを外す

すると,計算ディレクトリをWindows側からも見ることができるようになる.(図11)

Windows側からはファイルがフルコントロールになっているので,入力ファイルをWindows側のエディタで操作することもできる.(図12)これは,ちょっとした感動である.リンクを張っておくと,さらに操作性は向上する.もちろん,ファイルやディレクトリの名前変更や移動等も可能になっている.

図11 Windows から ubuntu のディレクトリが丸見え

図11 Windows から ubuntu のディレクトリが丸見え

図12 Windows のエディタで ubuntu の入力ファイルを直接編集することができる

図12 Windows のエディタで ubuntu の入力ファイルを直接編集することができる

同様に,Windowsファイルシステム内のファイルを読み込んで可視化することができる.先程はWSL内のV-Isioで可視化したが,今度はWindowsにインストールしたV-Isioで可視化してみる.V-Isioを立ち上げて,Windowsファイルシステム内のファイルを選択する.(図13)

すると,問題なくWindows側のV-Isioで可視化することができた.(図14)

当然だが図8と同様の可視化結果が得られている.

図13 Windows のV-Isioでubuntuの可視化データを直接読み込む

図13 Windows のV-Isioでubuntuの可視化データを直接読み込む

図14 2次元キャビティ問題における速度の絶対値のプロット図.Windows にインストールしたV-Isioで可視化.

図14 2次元キャビティ問題における速度の絶対値のプロット図.Windows にインストールしたV-Isioで可視化.

性能評価結果

これまでに,プリからポストまでの一連の流体解析を通じ,従来の仮想環境やエミュレータと同等以上の操作性を実感した.

とは言っても,計算性能に難があるようではやはり実用的には使い物にならないため,簡単にWSL上での計算時間や実行性能を調査することとした.

方法として,仮想環境(VMware)上にWSL環境と同一のバージョンのubuntuとFFV-Cを用意し,計算時間や計算性能を両者で比較した.ベンチマークテストとしては,2次元キャビティ問題を選んだ.総セル数は64×64×1=4096であり,圧力ポアソン方程式の反復法は点過緩和SORで,反復回数は最大20回とし,相対残差の閾値は10-4とした.なお,MPIプロセス数は1とし,スレッド数は4とした(図15中の(1)).

それぞれのプロファイルレポートを,図15と図16に示す.

Total execution time (マスタープロセスの経過時間)を両者で比較すると,後者が約3.80秒なのに対し前者は2.78秒で済んでおり,WSLでの計算が約1.37倍早いことが分かった(図15, 16中の(2)).総計算時間に対する各サブルーチンの割合を調べると,どちらの場合にも上位2つはProjection_Velocityと,Poisson_PSORであり(図15, 16中の(3)),これはどちらもフラクショナルステップの一部であり,ナビエ-ストークス方程式の演算部分にあたる.flop counts or byte counts 欄を見ると,どちらのサブルーチンもWSLのほうが1.1-1.2倍程度高い実効性能が出ていることが分かり,WSLは仮想環境と比較して独自のOSを立ち上げる必要がないためにCPUやメモリ等のリソースを有効活用されていることが示唆される.

その一方で,WSLでのFile_Outputはたった40回のコールで総計算時間の20パーセント弱を食っており,実効性能値を見てもVMwareでのそれと比較して1/10以下であることが分かった.他のファイルI/Oサブルーチンである History_outも同様の傾向を示しており,WSLでのファイルI/Oの性能がVMwareでのそれと比較して,極めて高コストであることを示している.

図15 プロファイルレポート(WSL)

図15 プロファイルレポート(WSL)

図16 プロファイルレポート(VMware)

図16 プロファイルレポート(VMware)

まとめ

Windows10上の ubuntu互換環境(Windows Subsystem for Linux, WSL)にオープンソースの流体ソフトウェアFFV-Cをインストールし,流体解析の一連の流れを通じてWSLの使い勝手を検証するともに,流体ソフトウェアの簡単な性能評価を行い,WSLにおける流体解析の可能性を調査した.

FFV-Cのインストールでは,必要とされるライブラリやソフトウェアがほとんどubuntuネイティブ環境下と同等であることが分かった.また,プリからポストまでの一連の流体解析も,従来のubuntuネイティブ環境下と全く同様に実行できることが確認できた.

もっとも,記事執筆時(2016/9/9)において WSLはベータ版であり,64ビット版限定での提供やシェルの日本語が一部文字化けする等,多少機能に不十分な点も見受けられた.しかし, Microsoft社の本気度(開発項目の優先順位などは[5]参照),導入の容易さ等から考えると,近い将来に従来の仮想環境やエミュレータを脅かす,ないしは置き換える存在になり得るものと考えている.

2次元キャビティ問題を対象にWSLの計算性能を調査した.その結果,WSLでの計算が従来の仮想環境での計算と比較して約1.37倍早いことが分かった.これは,WSLが仮想環境と比較して独自のOSを立ち上げる必要がないため,CPUやメモリ等のリソースが有効活用されていることを示唆している.一方で,WSLのファイルI/Oの実効性能値がVMwareでのそれと比較して1/10程度に留まり,WSLのファイルI/Oが極めて高コストであることが分かった.WSLでの計算では,以上の特性を留意しておく必要があると考えられる.今後は,プロセス/スレッド数依存性の調査や,他のオープンソース流体ソフトウェア(FrontFlow/BlueやOpenFOAM等)を対象とした同様の検証を行うことも検討したい.

近年,PCの性能向上や計算手法・格子生成技術の発達等により,ある程度大規模な計算がローカルPCでも可能になってきた.一般に業務用PCはWindowsマシンであるが,プリからポストまでの流体解析をほぼWindows操作することが可能なubuntu互換環境WSLを用いることで,従来の流体解析がさらに身近になったのではないかと考えられる.

 

参考文献およびウェブサイト

[1] https://msdn.microsoft.com/en-us/commandline/wsl/about

[2] http://www.vmware.com/jp.html

[3] https://www.virtualbox.org/

[4] https://www.cygwin.com/

[5] https://msdn.microsoft.com/en-us/commandline/wsl/faq

[6] http://avr-aics-riken.github.io/ffvc_package/

[7] https://github.com/avr-aics-riken/PMlib

[8] http://pc.watch.impress.co.jp/docs/column/nishikawa/1017333.html

[9] http://www.atmarkit.co.jp/ait/articles/1608/08/news039.html

[10] http://rcmdnk.github.io/blog/2016/06/05/computer-windows-ubuntu-bash/

[11] http://cabonera.hateblo.jp/entry/2016/08/04/003300

[12]http://answers.microsoft.com/ja-jp/windows/forum/windows_10-performance/anniversary-update/66e03b07-4845-4a66-be07-0feda12bd34e

[13] https://www.microsoft.com/ja-jp/software-download/windows10

[14] https://www.open-mpi.org/

[15] User Guide of FFV-C Frontflow/violet Cartesian

[16] http://mobaxterm.mobatek.net/

[17] http://avr-aics-riken.github.io/V-Isio

『分野4次世代ものづくり』シンポジウム最終成果報告会レポート

 去る3月23、24日に行われた文部科学省HPCI戦略プログラム『分野4次世代ものづくり』シンポジウム〔最終成果報告会〕のレポートをお届けします。
6年間に及んだ研究開発の成果に関する詳細な報告を求める方は今後発表される別の資料も参照してください。計算工学ナビ・ニュースレター(ダウンロードページ)でも加藤千幸教授による総括を読むことができます。ここでは、シンポジウムで感じることができた『次世代ものづくり』の今と未来の可能性について概要をお伝えします。

第1日

成果の概要とその先への展望

sympotop まず、分野4統括責任者の加藤千幸教授(東京大学生産技術研究所 革新的シミュレーション研究センター長)より、プロジェクトの所定の目標を達成したことが宣言されました。研究課題ごとの達成状況を見ると「大幅に達成」と評価された研究項目がいくつかあります。これはものづくりのイノベーションにつながる分野4の特徴的な成果に対する評価であると筆者は捉えました。そのあと、この2日間のシンポジウムで発表される各課題の概要が示され、最後に「これはスターティングポイントであり、今後につなげることが重要。このシンポジウムでは近未来のものづくりを実現するHPCについて議論したい」と、会の位置づけを明らかにして締めくくり、各研究課題の担当者による成果発表へと移りました。

課題1:輸送機器・流体機器の流体制御による革新的高効率化・低騒音化に関する研究開発

sympo_k1f 課題責任者の藤井孝藏教授(東京理科大学、JAXA宇宙科学研究所)は、この研究を通じて示したかったのは「計算機の中で新しいアイデアを試す」ことである、という話から発表を始めました。翼の形状の工夫によって性能を向上させることが限界に達してしまった現状を、プラズマアクチュエータという動的に機能する小さなデバイスによって打開することがこの課題の目的ですが、それに必要な知見を計算機(HPC)の活用によって得ることが重要ということでしょう。
藤井教授のチームは京コンピュータを活用して、低・中レイノルズ数領域における現象理解を深め、デバイスの効果的な利用方法を見つけました。この知見を設計の指針として提供できたことが、最大の成果であると藤井教授は位置づけています。そして、この成果のものづくりへの適用が企業との共同研究を通じて始まっています。
次の講演者、株式会社東芝の田中元史参事は、同社のメガワット級発電用風車にアクティブ流体制御技術を適用する実証研究について発表しました。変動する自然風のなかでも流れを制御できる見通しが得られつつあること、そして、すでに鹿児島の2MW風車を用いて約1年間のフィールド実験を進めてきたことが公開され、プレゼンテーションの最後には、巨大な風車をドローンで撮影した雄大な映像が紹介され、会場から大きな拍手が起こりました。
共同研究事例の2社目として登壇したのは、マツダ株式会社の清水圭吾氏です。プラズマアクチュエータをクルマの車体に装着し、デザイン自由度を損なわずに空力性能を高めるという野心的な研究です。具体的には、デザイン上重要な車体後部の丸みを残したまま、プラズマアクチュエータによって気流の剥離を促してドラッグを軽減するというアイデアの有効性を、京を使った大規模流体計算によって確認しました。冒頭で藤井教授が述べた、「計算機の中で新しいアイデアを試す」ことの威力が証明された事例であると感じました。
本課題最後のプレゼンテーションは、産業技術総合研究所瀬川武彦主任研究員による、実験研究者の視点によるプラズマ気流制御技術に関する現状の報告でした。瀬川氏らは、一般的なシート型プラズマアクチュエータよりも設置が容易で空力性能に影響を与えにくい、ひも型プラズマアクチュエータを開発しています。分野4の成果として示された剥離制御に関するデータや知見が、実験データとの比較対象としてきわめて有効であること、無限にあるパラメータの最適化に役立つことが発表されました。分野4の取り組みが、志を同じくする他の研究グループにも利用されている事例として興味深く見ることができました。

課題2:次世代半導体集積素子におけるカーボン系ナノ構造プロセスシミュレーションに関する研究開発

sympo_k2 シリコンカーバイド(SiC)やグラフェンといった新しい材料を採用したデバイスの開発に、従来のシリコンベースの技術を適用するのは難しいため、今まさに新しい知見が求められていると、課題責任者の大野隆央氏(物質・材料研究機構NIMS特命研究員)は研究開発の背景を説明しました。大野氏らのグループは、第一原理分子動力学計算プログラムPHASE/0を京に対して最適化して、高精度長時間の解析を繰り返すことで、グラフェン成膜プロセスとSiC酸化プロセスに関する世界初と言っていい画期的な知見を得ることに成功しました。それぞれの成果に関する3名の講演者による発表の概略を順に紹介します。
NIMS特別研究員の山崎隆浩氏は、高品質なグラフェンの生産手法=レシピを確立するために必要な、SiC表面上での個別原子レベルの成長機構を解明するために、PHASE/0と京を用いて原子数900〜2000規模のシミュレーションを実施し、その結果、炭素の鎖が成長・移動して2次元的構造を形成することを示しました。プレゼンテーションでは、精細なアニメーションが上映され、動画ベースの観察が現象の理解に効果的であることが強調されました。
続いて、NTT物性科学基礎研究所の高村真琴氏は、シミュレーションと実験の両輪で研究を進めることの重要性を、SiC上グラフェンの成長メカニズムを例に指摘しました。特に成長の初期過程でどのようなメカニズムが作用しているかを知る上で、今回得られた知見は重要だったとし、今後はシミュレーションによって実験の指針を得ることができるのではないかという期待を明らかにしました。
株式会社東芝研究開発センター研究主幹の清水達雄氏からは「SiCパワーデバイス特性向上を目指したシミュレーション研究」と題して、市場の拡大が進む次世代パワーデバイスの重要性と開発における課題についてわかりやすい解説がありました。現在進められているSiCベースのパワーデバイス(縦型MOSFET)の開発にはまだ技術課題がいくつか残っていますが、大規模シミュレーションによって改善のヒントが得られることが期待されており、今回NIMSと東芝の共同研究から得られた成果もそのひとつとなることが望まれています。

課題3:乱流の直接計算に基づく次世代流体設計システムの研究開発

このセッションでは、自動車、船舶、燃焼器といった複雑な工業製品の設計に京クラスのHPCを導入することでどのようなイノベーションが起こるかが、多くの具体例とともに示されました。特に自動車会社10社、サプライヤー4社、大学研究7機関が参画するコンソーシアムによる、自動車関連の実証研究の成果が多く発表されたので、まずそこから紹介します。
sympo_k3 まず、神戸大学の坪倉誠教授からは、0.5mmを解像する極めて大規模な空力解析により、2%以内の誤差で自動車の風洞実験を代替可能という明確な成果が発表されました。風洞では実現が難しいレーンチェンジによる蛇行や、風がクルマに与える影響などを加味した評価が行える点も、この技術の革新性です。
スズキ株式会社の飯田桂一郎氏からは、ADVENTUREclusterやFrontFLow/blue-Acousticsといったソフトウエアを使って連成解析を行うことで、車室内の騒音を良好な精度で予測できることが示されました。設計上「一番知りたい」とされる、4KHz以下の周波数特性が予測できたことから、今後の実用化が期待されます。
株式会社本田技術研究所主任研究員の寺村実氏からは、「Hondaでの京を活用したアクティビティ」として、ミラー騒音解析の事例が紹介されました。ミラー単体に対して4.2億セルというメッシュ解像度を適用することで、実験値に近い流れ現象の解析が可能であることが示されました。
富士重工株式会社の小林竜也氏からは、「車両運動性能向上を目指した非定常空力シミュレーション」というテーマで、従来は困難だった、高速走行時に揚力の低周波振動を生み出す流れの評価を、スーパーコンピュータによる非定常シミュレーションで行う事例が紹介されました。
続いて、一般財団法人日本造船技術センターの西川達雄課長より、次世代数値曳航水槽の現状に関する報告がありました。すでに、従来型の曳航試験に対して誤差1%以下という極めて高精度なCFD計算が可能になっています。今後は、実用化を加速するための研究開発を進めていくと述べました。
この研究開発課題最後の発表は、京都大学・武藤昌也特定助教による燃焼器設計支援のための燃焼シミュレーションに関する発表でした。ガスタービン燃焼器、環状燃焼器、発電用のマルチバーナ微粉炭ボイラなど、様々な実機を対象に大規模計算を実施した結果得られた知見が紹介され、産業界での応用に向けた手応えが共有されました。

第2日

課題4:多目的設計探査による設計手法の革新に関する研究開発

多目的設計探査とは、多目的進化計算の結果に基づいて多ケースのシミュレーションを実行し、設計における複雑なトレードオフの関係を明らかにする手法です。
JAXAの大山聖准教授を中心とするチームは、大規模な設計問題にこの手法を適用するため、アルゴリズムとソフトウエアの開発を進めてきました。京コンピュータを使うことで初めて可能となった成果が次々とあがっています。
sympo_k4 東京理科大学の立川智章講師からは「京で可能になったロケット射点形状の空力音響多目的設計探査」と題して、ロケットの打ち上げ時に超音速ジェットから発生する強い音響波の影響を射点(発射台)の形を適切に設計することで低減する試みとその成果について発表されました。射点の形状をパラメータ化し、多目的進化計算に基づいて異なる形状を生成して、それぞれについて流体ミューレションを実行することで最適解を探索します。京の6500ノード(約7%)を用いて約2週間の計算を実行した結果、音響強度最小化、表面圧力最小化および形状簡略化の間のトレードオフを明らかにしました。そのプロセスをふんだんなビジュアルを使って説明するプレゼンテーションは興味深く見応えがありました。
続いての発表は企業における製品開発への適用です。横浜ゴム株式会社の小石正隆理事・研究室長によるプレゼンテーションは、多目的設計探査の大きな可能性をわかりやすく伝えるものでした。第1世代のフィンタイヤ(タイヤ側面にフィン状突起を設け空気抵抗を低減した自動車タイヤ)を改良して、空気抵抗だけでなくリフト(揚力)も低減する第2世代フィンタイヤを開発するために多目的設計探査の手法とFrontFlow/redが使用されました。事前計算を含め100万ノード時間を超えるパラメトリックスタディを実施した結果、抵抗とリフトを同時に低減できる革新的な知見が得られ、それを裏付ける実車試験結果も得ることができました。プロトタイプは2015年の東京モーターショーにおいて公開され、注目を浴びました。
マツダ株式会社の小平剛央アシスタントマネージャーからは、「スパコン京を用いた複数の車体構造の同時多目的設計最適化」というテーマで先進的事例が紹介されました。同社では複数の車種の基本構造となるコモンアーキテクチャ構想を進めています。従来の最適化技術は単一車種にしか対応していないため、同社とJAXAの共同研究によって複数車種同時最適化という課題に取り組みました。京の525万ノード時間を使って16,320回の衝突シミュレーションを行う過程で、「予測はしていたが本当にあるのかわからなかった解」が見つかったといいます。将来は新しい商品の創造にも適用したいと、この技術の可能性に対する高い評価が印象的でした。

課題5:原子力施設等の大型プラントの次世代耐震シミュレーションに関する研究開発

sympo_k5 課題責任者である日本原子力研究開発機構(JAEA)の中島憲宏副センター長から、耐震シミュレーション技術の京での利用とその高度化状況、具体的なプラントや機器を対象とした数値実験などの成果について報告がありました。JAEA内の課題に対する取り組みとしては、高温工学試験研究炉(HTTR)に対する実証計算の事例が紹介されました。約30万点の部品からなる大規模な組立構造物の解析を可能とし、HTTRの観測値と比較することでシミュレーションの精度が検証されました。
産業界との共同研究については、千代田化工建設株式会社の松川圭輔主席技師長から「石油・化学プラントにおける次世代耐震シミュレーションの活用」、株式会社荏原製作所の杉山道子構造・振動解析グループ長から「一般産業機械における次世代耐震シミュレーションの活用」 というテーマで、京によって可能になった耐震シミュレーションの世界が紹介されました。

計算科学技術推進体制の構築

sympo_h 5つの研究開発課題に関する発表のあと行われたのが、計算科学技術推進体制の構築に関するセッションです。責任者の畑田敏夫特任教授(東京大学生産技術研究所)は、この事業の取り組みを「(プロジェクトの成果を)使いやすくする」「知ってもらう」「活用できる人材を育てる」「実際に使ってもらう」という4つにカテゴリーに分類して説明しました。産業界への貢献を重視する分野4の特徴的な事業と言えるでしょう。
成果を使いやすくすることを目的とした具体的な取り組みとして、次に理化学研究所 計算科学研究機構の小野謙二チームリーダーよりHPCプラットフォーム(HPC/PF)の開発と成果に関する発表が行われました。HPC/PFは最先端の大規模並列シミュレーション技術を研究・設計の現場で使える「道具」にする基盤ソフトウエア群の名称です。一般的なPCクラスタから京コンピュータまでの多様な環境で運用でき、ソフトウエアだけでなく解析事例のデータベースも提供するこのプラットフォームは、すでにハンズオンセミナー等を通じて体験可能な状態です。
一連の発表の最後に、川鍋友宏特任研究員(東京大学生産技術研究所)からアウトリーチ活動の成果について報告がありました。アウトリーチ活動の中心的な役割を果たしているのが、いま皆さんがご覧になっているWebサイト『計算工学ナビ』です。2015年3月までの訪問者は43,000人、印刷版ニュースレターの発行部数は10,000部となりました。また、解析事例データベースには200件を超えるデータが登録済みで、簡単なユーザー登録によって全件を検索・閲覧することが可能です。今後も情報の拡充を図っていきます。

パネルディスカッション「今後の実用化に向けて」

2日間のシンポジウムの最後に、下記のパネリストを迎えて、今後の実用化をテーマにパネルディスカッションが行われました。後半には会場からも多くの意見が出され、大変活発な議論となったことが印象的でした。

sympopanel司会:
加藤 千幸(東京大学生産技術研究所 教授)
パネリスト:
藤井 孝藏(東京理科大学 教授)
大山 聖(宇宙航空研究開発機構 宇宙科学研究所 准教授)
西川 達雄(一般財団法人日本造船技術センター 課長)
農沢 隆秀(マツダ株式会社 技術研究所 技監)
福島 伸(株式会社東芝 研究開発センター 首席技監)

この記事では、各パネリストの発言について個別に紹介することはしませんが、加藤千幸教授による冒頭言について簡単に紹介したいと思います。
京を活用する高性能なアプリケーションソフトを開発し、企業との共同開発等を通じてその有効性を実証するという当初の目的は十分果たしたと加藤教授は述べました。次に「本当に京でしかできない革新的な成果はできたのか?」と自問し、これについても「できた」と答えます。シンポジウムで発表された印象的な成果として、JAXAと東芝による風車へのプラズマアクチュエータの適用、NIMSと東芝によるSiC酸化プロセスのメカニズム発見、日本造船技術センターの数値曳航水槽、1ペタフロップス級のスパコンがあれば大規模設計探査が可能であるとするマツダの検証結果、千代田化工による耐震裕度の解析技術などに触れました。しかし、一方で「スッキリしないものが残っている」とも発言し、会場を刺激しました。開発したアプリがすぐさま実用化できるレベルにないことが、その理由です。ポスト京に進む上で、今後どのような体制が必要か、もういちどよく考えることが重要だと呼びかけました。

レポート:計算工学ナビ編集部

ABINIT-MP Openシリーズ

abopen00
 
2016年12月版に関する新しいページがあります

はじめに

logoABINIT-MPは、フラグメント分子軌道(FMO)計算を高速に行えるソフトウェアです[1]。専用GUIのBioStation Viewerとの連携により、入力データの作成~計算結果の解析が容易に行えます。4体フラグメント展開(FMO4)による2次摂動計算も可能です。また、部分構造最適化や分子動力学の機能もあります。FMOエネルギー計算では、小規模のサーバから超並列機の「京」まで対応しています(Flat MPIとOpenMP/MPI混成)。
 

特徴

ABINIT-MPは使いやすいFMOプログラムで、4体フラグメント展開までが可能です。研究室単位のLinux/Intel系サーバに標準搭載されているMPI環境で動作しますし、特別な設定も必要ありません。また、煩雑で注意深さを要するフラグメント分割を伴う入力データの作成は、随伴GUIのBioStation Viewer(Windowsで動作)を使うなどすれば容易に作成できます。また、フラグメント間相互作用エネルギー(IFIE)などの計算結果は膨大となりプリントからの理解はしばしば困難ですが、Viewerを使うと可視的・直観的に対象系の相互作用の様態を把握できます。

主な機能

・エネルギー
→ FMO4: HF, MP2 (CD)
→ FMO2: HF~CCSD(T)
→ FMO2: CIS, CIS(D)系
・エネルギー微分
→ FMO4: HF, MP2
→ FMO2: MP2構造最適化, MD

続きを表示

・その他機能
→ SCIFIE, PB, CAFI, FILM, α(ω)
→ GUI(BioStation Viewer)
→ MPI, OpenMP/MPI混成
→ 最深部はBLAS処理

■BioStationViewer / GUI
abopen01
 

開発の経緯

ABINIT-MPプログラムは、東京大学生産技術研究所を拠点とする「戦略的基盤ソフトウェアの開発」、「革新的シミュレーションソフトウェアの開発」、「HPCI戦略分野4 次世代ものづくり」の一連のプロジェクト(代表:東京大学 加藤千幸教授)、さらにJST-CREST「シミュレーション技術の革新と実用化基盤の構築」(代表:神戸大学田中成典教授)、科学研究費補助金:特定領域研究「実在系の分子理論」(代表:京都大学 榊茂好教授)、立教大学SFR(担当:立教大学望月祐志教授)などの支援を得て、10年以上に渡って開発が進められてきました。Intel IA64系バイナリは、Ver.7が東京大学生産技術研究所の革新的シミュレーション研究センターで、また「京」向けのVer.6+が理化学研究所計算科学研究機構で利用可能となっています。

今後は、東京大学工学研究科を代表拠点とする「フラッグシップ2020 重点課題6」(代表:東京大学 吉村忍教授)の中でものづくり系への応用を指向して改良とリリースが行われていく予定です。後述のように版管理の枠組みも変わってOpenシリーズに移行します。

ABINIT-MPの開発とその先導的な実証応用には多数の機関から多くの方が関わってきましたが、物理的な活動の場所としては東京大学生産技術研究所の一角に置かせていただいてきました。こうした経緯もあり、ABINIT-MPの新しいポータルHPを「計算工学ナビ」の中に設けていただくことになりました。

開発者(所属)

望月祐志*(立教大学 理学部), 中野達也(国立医薬品食品衛生研究所 医薬安全科学部), 坂倉耕太(NEC), 山本純一(NEC), 沖山佳生(元 東京大学 生産技術研究所), 山下勝美(元 NECソフト), 村瀬匡(元 NECソフト), 石川岳志(長崎大学 医歯薬学総合研究科), 古明地勇人(産業技術総合研究所 バイオシステム部門), 加藤雄司(元 立教大学 理学部), 渡辺尚貴(みずほ情報総研), 塚本貴志(みずほ情報総研), 森寛敏(お茶の水女子大学大学院 人間文化創成科学研究科), 田中成典(神戸大学大学院 システム情報学研究科), 加藤昭史(みずほ情報総研), 福澤薫(日本大学 松戸歯学部), 渡邉千鶴(元 東京大学 生産技術研究所)
(*取り纏め責任者)

応用分野

ABINIT-MPのFMO計算は、開発当初から生体分子関係、特にタンパク質とリガンド(薬品分子)の複合系に対して主に用いられてきました。これは、計算で得られるフラグメント間相互作用エネルギーがアミノ酸残基間、あるいはリガンド-アミノ酸残基間の相互作用の状態を理解するのに好適なためです[1,2]

続きを表示

しかし、FMO計算は生体系に限られるだけでなく、水和凝集系や一般の高分子、あるいは固体なども扱える潜在力を持っています。ABINIT-MPの応用は、今後はこうした一般の化学工学や材料科学、あるいは応用物理関係の分野へも広がっていくことを期待していますし、そのための整備とエビデンスの集積も進めています。

abopen03
 

FMO創薬コンソーシアム

2014年11月から「FMO創薬コンソーシアム」(代表:日本大学 福澤薫助教)が産官学で組織され、「京」を計算資源としてABINIT-MPによるFMO計算に基づくタンパク質・リガンドの相互作用解析が進められています。重要な創薬ターゲットが設定されており、当該領域の共有基盤となる知見(特にIFIEのデータセット)の蓄積が期待されます。ABINIT-MPは改良が続けられていきますが、このコンソーシアムは実践的な利用者コミュニティとして今後も重要な役割を果たしていくことになります。

Openシリーズと今後のリリース

ABINIT-MPには、現在からバイナリで利用可能なVer.7(東京大学生産技術研究所からのインテル向け)とVer.6+(理化学研究所の「京」向け)とは別に、主に開発経緯的な事由から「ローカル版」が存在しています[1]。そちらでは、励起エネルギーや動的分極率の算定、さらに結合クラスター展開による高精度エネルギー計算などが利用できます。こうした機能に関心を持たれる方も居られますし、手持ちのマシン環境によっては再コンパイルやチューニングのためにソースを所望される場合もあります。こうした状況を改善すべく、産官学を交えたコンソーシアム的な組織でABINIT-MPのソース共有を行い、継続的なコード開発・改良と保守を図っていく活動の中でリリースされていくのがOpenシリーズになります。

続きを表示

Ver.7に準拠し、モデル内殻ポテンシャル(MCP)機能を追加した最初のOpen Ver.1については準備が整いました(2016年3月31日)。また、「フラッグシップ2020 重点課題6」のプロジェクト活動の中で2016年度に整備するOpen Ver.2では、密度汎関数(DFT)のモジュールなどを分子科学研究所の石村和也研究員のSMASHから、また2電子積分の恒等分解(RI)のモジュール(C言語で記述)を長崎大学の石川岳志准教授のPAICSから移植する他、「ローカル版」のいくつかの機能を含める予定です。また、ホットスポットの洗い直しやメモリ割り当ての見直しなども行い、パフォーマンスの総合的な向上も図ります。Open Ver.3以降も機能充実を順次進めていきます。また、BioStation Viewerについても随時Openシリーズへの対応を図ります。

ソースの共有とは別にOpenシリーズでも従来のバイナリでの提供も続けるつもりですが、これまでのIntel IA64系と富士通系だけでなく、NECのベクトル系(SX-ACE)、さらにパソコン用にWindows 64ビット系を提供しようと考えています。準備ができましたら、このページで順次ご案内します。

開発系コンソーシアム

ABINIT-MPのOpenシリーズの開発や保守にソースレベルでコミットしていただくための産官学枠組みです。バグ情報と対策、新規開発の機能のシェアなど意図していますが、参画される企業様が商用に独自の高速化や改良を図ることは基本的に可とする方針です。現初期段階では、個別にご参画をお願い・確認させていただいて立上げようとしているところですが、規模的には産官学で20拠点程になります。今後、このコンソーシアムについても情報を更新していく予定です。

コンタクト

ABINIT-MPのOpenシリーズ、あるいは開発系コンソーシアムにご関心のある方は、立教大学の望月(fullmoon -at- rikkyo.ac.jp)、中村(mnakamura -at- rikkyo.ac.jp)にメールでご連絡いただければと思います。

[1] “Electron-correlated fragment-molecular-orbital calculations for biomolecular and nano systems”, S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys. 16 (2014) 10310-10344.
[2] “The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems”, edited by D. G. Fedorov, K. Kitaura, CRC Press, Boca Raton, Florida, 2009 ISBN 978-1-4200-7848-0.

RaspberryPiクラスタ製作記 第4回「BareMetal並列計算」

最終回 BareMetal環境で並列計算をしてみよう

manga4

前回はOSの実行にかかるコストを削減するため、BareMetal環境で開発を行いましたが、実行時間はLinux上での結果に負けてしまいました。悔しいです。何とかしてリベンジできないでしょうか。

前回を思い返してみると、Linux上でmpichを用いてLU分解を並列化した結果は、ノード数を増やせば増やすほど通信コストの増加により実行時間は長くなっていました。 ということは、BareMetal環境を活かして通信コストを削減できれば、並列計算時の実行時間は勝てるかもしれません。

そこで今回は、このアイデアを実行し、リベンジを図りたいと思います。また、本連載は今回が最後なので、 記事の終わりに今までの振り返りと要点をまとめます。

通信コスト削減による高速化

冒頭のマンガで述べられているように、BareMetal環境は何もないので、 Ethernetを用いたネットワーク通信を行うには、 LANコントローラのドライバを書き、TCP/IPプロトコルスタックを実装しなければなりません。 この作業を行うとそれだけで1冊本が書けてしまうので、今回はもっと簡単な通信を用います(この後マンガの主人公は、チップの仕様書とUSBの規格書、TCP/IPプロトコルの解説書とにらめっこする日々を過ごすことになるでしょう)。

一般にBareMetal環境やマイコンなどで用いられる通信方法には、以下の3つがあります。

・UART
・SPI
・I2C

これら全てに共通するのは、通信仕様がUSBなどに比べ非常に簡単であり扱いやすく、通信にかかるコストが非常に少ない点です。 この点を活かせば、BareMetal環境上で通信コストを削減した並列計算が行えるはずです。

Raspberry PiのSoCに実装されている中で、最も速度の出る通信方式はSPI(Serial Peripheral Interface)です。 SPIは親(Master)と子(Slave)という役割に分かれ、4本の線を用いて行うシリアル通信です。

SPIconnect

通信は、マスターからスレーブセレクト線(SS)で選択されたスレーブのみが、クロック線(SCLK)から供給されるクロックに合わせて、1bitずつデータの送受信を行います。

SPIwave

SPIの通信クロック(SCLK)は、ペリフェラルマニュアルの156ページより、

SCLK = Core Clock(default: 250MHz) / CDIV(=2^n)  # ただしn=0の場合はCDIV=65536

で計算されるので、n=1の場合に最大の125MHz、つまり理論的には最大125Mbpsで通信が行えます。
しかし、私がRaspberry Pi同士のSPI通信を試したところ、Slave側がうまく動作せず、通信できませんでした。問題解決のため、配線、マスタ側からの出力、Slave側のGPIO設定、SPI Slaveレジスタ設定などの確認を行いましたが問題はなく、加えて受信処理プログラムのミスを確認するためSPI Slaveテストレジスタを用いたテスト通信を行いましたが、こちらも問題なく動作しました。つまり、Slave側でIOポートを通した通信を受信することができていないようですが、今のところ私にはなぜうまくいかないのかわかりません。この問題については引き続き調査を続けますが、記事の締め切りがあるため今回は別の方法を用いることにしました。

2番目に高速なのはUART(Universal Asynchronous Receiver Transmitter)です。 これは通常9600 baudや115200 baud程度の低いボーレート(1秒間に行われる変復調の回数)で、パソコンやマイコン間で文字列を送受信するために使われていますが、設定を変更することでボーレートをMbaudまで上げた高速通信が可能です。

代わりにSPIと違い、基本は1対1での通信しか行えません。 この点は非常に残念ですが、ほかに高速で簡単に扱えるものが見つからないためBareMetal環境の通信はUARTを用いて行います。

Ethernet&TCP/IP vs UART

BareMetal環境の通信にはUARTを用いることにしましたが、 理論的にはEthernet&TCP/IPと比べてどれだけ速くなるか気になります。

そこで、ここではEthernet&TCP/IPとUART&独自プロトコルの速度を調べて比較します。

Ethernet&TCP/IP

ご存知のように、ネットワーク通信は送受信するデータを細かく分けたパケットという単位で行います。 このパケットは入れ子構造になっており、例えばTCP通信のパケットは、Ethernet, IP, TCPパケットに包まれた3重構造になっています。

各パケットには宛先などの情報が記されたヘッダがついており、Ethernet、IP、TCPパケットのヘッダを合わせると、長さは66byteになります(オプションなしの場合)。 一番外殻であるEthernetパケットの長さは最大で1526byteなので、 Ethernetを用いた通信は、データ1460byteにつき66byteのオーバーヘッドが発生します。 このオーバーヘッドを考慮して最大通信速度を考えると、100Mbps * 0.96 = 96Mbpsとなります。

UART

Raspberry PiのUARTのボーレートは、ペリフェラルマニュアルの183ページより、次の計算式で求められます。

baudrate = F_UARTCLK(default: 3MHz) / (16 * BAUDDIV)

この計算式のうち、F_UARTCLKとBAUDDIVは変更可能であり、前者は起動スクリプトを修正、後者はレジスタの設定を行うことで変更できます。

F_UARTCLKは通常3MHzとなっていますが、SDカードの第1パーティションにあるconfig.txtに次の一文を追記すると変更できます。

init_uart_clock=416000000

ここに設定する値は、Raspberry PiのSoCが搭載しているUARTペリフェラル(PL011)マニュアルによると、Core Clock(250Mhz)* 5/3 までの値が推奨されているので、F_UARTCLKは416MHzが最大となります。

BAUDDIVはIBRDとFBRDレジスタで設定でき、 どちらも0にした場合にBAUDDIV=1となり、 最大ボーレートは以下の計算式より26Mbaudとなります。

baudrate = 416(MHz) / (16 * 1) = 26(Mbaud)

UARTは1度に1bitの変復調を行うので、通信速度の最大は26Mbpsとなります。

実際に設定してオシロスコープで見てみると、確かに26MHzの波形が出ていることが確認できます。

UART_freq_26MHz_recvtest

UARTの通信は8bitずつ行われ、データの開始と終了を判定するためにスタートビットとストップビットが最小で1bitずつ付きます。

UARTsignal

つまり、8bitあたり2bitのオーバーヘッドが発生するので、これを考慮した最大通信速度は 26Mbps * 0.8 = 20.8Mbpsとなります。

Ethernet&TCP/IPとUARTの比較

以上で述べた、EthernetとUART通信のオーバーヘッドを含めた最大速度を以下に示します。

通信方式 最大通信速度(Mbps)
Ethernet&TCP/IP 96
UART 20.8

この結果はEthernetの圧倒的勝利のように見えます。しかし、現在考慮しているEthernet&TCP/IPのオーバーヘッドはトランスポート層までであり、セッション層以降のオーバーヘッドについては考慮していません。 つまり、セッション層以上にまだまだ多くのオーバーヘッドが存在するはずです。 対するUARTは1対1通信のため送信先選択などを行う必要がなく、オーバーヘッドは少なくなるはずです。 この差を考慮すれば、もしかするとUARTを用いた並列計算のほうが速くなるかもしれません。

まだ希望はあるはずです。 やって確かめてみましょう。

ソフトウェアの実装

UARTは1対1通信なので、前回作成した並列計算プログラムをノード数2に限定し、行列Aを2つに分けた前半を処理するノード0と、後半を処理するノード1のコードを作ります。

MPI関数を用いて行っていた枢軸ベクトルや計算後のLU行列の転送は、以下のUART通信を用いた関数に置き換えて行います。

void Serial_send(uint8_t *buf, int len) 第1引数bufで指定されたアドレスから第2引数len(byte)データを送信する
void Serial_receive(uint8_t *buf, int len) 第1引数bufで指定されたアドレスから第2引数len(byte)データを受信する
void uart0_putc(int c) 第1引数cで指定された値の下位8bitを送信する
int uart0_putc() 受信した8bitの値を返す

この変更を行ったコードをこちらで何度か動かしたところ、UARTで10byte程度以上の連続通信を行うと、データの整合がとれなくなる問題が発生しました。そのため、データの転送は4byte(double型の1つ分)ごとに 開始文字’s’と終了文字’e’を互いに送り合い、同期を取るようにします。

以上を行ったコードが以下の2つになります。

baremetalLU_node0.zip
baremetalLU_node1.zip

今回は通信にUARTを用いるため、UARTを用いた計算結果や実行時間の出力が行えません。そこで出力はSPI通信で別のマイコンに転送し、そのマイコン経由でパソコンに表示します。

マイコンには開発の容易さからmbed LPC11U24を用います。マイコンのソースコードは以下からダウンロードして使用してください。

spi2serial_full

なお、mbed LPC11U24のボーレートは115200bpsに設定しています。

ハードウェアの接続

データの通信を行うUARTと、結果の出力を行うmbed LPC11U24の接続を以下に示します。

raspi_uart

この図を参考に配線を行ってください。
なお、ログ出力用のマイコンはUSBケーブルを用いてパソコンと接続してください。

動作比較

プログラムのコンパイルを行い、出てきたバイナリをRaspberry Piにセットし、電源を入れるとプログラムが動作します。果たして、通信にUARTを用いた結果、通信コストの削減とそれにともなう処理の高速化はできたのでしょうか。以下に結果を示します。

exec time: 6278386

結果は6278ミリ秒。

Linux上で2ノード並列計算を行った際の結果は1600ミリ秒程度だったので、約4倍も遅くなってしまいました。通信速度で既に大きな差があった上に、データの同期を取るため8byteにつき2byteの大きなオーバーヘッドを作ってしまったことが敗因でしょう。

結局、私によるBareMetalの挑戦は残念な結果に終わりました。
もちろん、Linuxも同じRaspberry Piで動いているのですから、同じ設定をハードウェアに行えば、BareMetal環境でも同等かそれ以上の性能が得られるでしょう。しかし、そのためにはより多くの時間と労力が必要になり、難しいです。

片桐さんの著書「スパコンプログラミング入門 並列処理とMPIの学習」の35ページには、次のように書かれています。

「実装する処理によっては,スパコンのメーカー製やフリーソフトウェアによる数値計算ライブラリが利用できることがあります.この場合は,自らプログラミングするのではなく,提供されている数値計算ライブラリを利用する方が,性能が格段に良いことが多いです.〜中略〜最適化されたBLASの性能は,個人で開発した行列と行列の積に対して,10倍以上高速であることも珍しくありません」

私の挑戦が失敗したことからわかるように、既存のOSや計算ライブラリは既に高度な最適化を施されており、高速に動作します。そのため、計算処理の高速化は既存のOSや計算ライブラリを基にチューニングを行うほうが、手軽に期待する性能が得られるでしょう。しかし、そのチューニングのためには、ハードウェアやCPUアーキテクチャについての知識が必要となります。BareMetal開発はこの知識を得るために行い、そこで得た知識を既存のものに対し適応するという手法が良いのではと考えます。

連載のまとめ

本連載は今回が最後なので、今までのまとめとふりかえりを行います。

第0回

第0回では、

・そもそもスパコンとはなにか
・スパコンで計算される問題は本当にスパコンを用いる必要があるのか

を確認するため、現在の一般的なパソコンとスパコン「京」の性能を比べ、 実際に「京」で使われているPHASE0というプログラムを動かして、 実行時間が指数関数的に増えていく様子をグラフで確認しました。

現在はパソコンの性能が上がり、以前はスパコンを用いていた計算も、パソコンで行えるようになりました。 そのため、大学でコンピュータを学んでいる人の中でも、実際にスパコンに触れられる方は少なく、 スパコンとはなにか、何故必要なのかがわからない方もいるのではと感じています。 しかし、調べてみるとスパコンが必要となる計算は確かに存在し、 それらは私たちの生活に深く関わっています。 それを少しでも感じていただけたなら幸いです。

第1回

第1回では、

・クラスタとはなにか
・自作スパコンの作り方

について説明しました。

現代のスパコンは、複数のマシンをつなげてクラスタを作り、高い性能を得ています。 この「クラスタマシンを作る」というと、実用重視で高性能高価格なマシンを繋いだものを想像し、個人では手を出せないと考えてしまうかもしれません。 しかし、勉強のためと割り切って考えれば性能は必要なく、今回のように安価なRaspberry Piを使って作ることも可能です。

また、今回は16台でクラスタマシンの制作を行いましたが、 並列計算の学習を行うだけであれば4台もあれば十分でしょう。 そうすると、コストは1/4の3万円程度になるので、手を出しやすいかと思います。

第2回

第2回では、第1回で制作したRaspberry Piスパコンを用いて、

PAHSE/0の動作確認
・Linpackを用いた性能評価
・スパコン「京」との性能&コスト比較

を行いました。

PAHSE/0 と Linpack、どちらもノード数を増やすごとに処理速度は上がりましたが、その伸びは非常に悪くなっていきました。 ノード数(プロセッサ)と処理速度の関係には「アムダールの法則」という法則があり、 この法則よりノード数を増やせば増やすほど並列処理に含まれる並列化できない部分がネックになり、 速度の向上の伸びは非常に悪くなっていくことが示されます。

第3回

これまでは既存のソフトウェアを用いて性能評価などを行いましたが、 第3回では実際に自分でMPIを用いて並列処理を書いてみました。

これまでのPHASE/0とLinpackを実行した結果は、プログラムの並列化を行えば(伸びが悪くなることはあっても)必ず高速化されていました。 しかし、自分で並列化したコードはノード数を増やすごとに実行時間が増えるという結果になり、 並列化によりきちんと高速化されるプログラムを作るのは難しいことがわかりました。

今回用いた初代Raspberry PiはCPU性能が低く、並列化のためのコストが結果に大きく響きます。 逆に言えば、正しく効率的な並列化を行った場合のみ性能が向上する素直なマシンであり、 分散メモリ型の並列プログラミングの教材として最適なマシンと言えます。

一方、2015年2月3日に発売された新型のRaspberry Pi 2は、CPUが ARM11 700MHz シングルコア から Cortex-A7 900MHz クアッドコア(4コア)に変更され、性能が大幅に向上しました。発売直後、個人的に購入した1台のRaspberry Pi 2上でコアを4つ用いてLinpackを実行してみたところ、1台で旧型(Raspberry Pi Model B+)6台分の性能が得られました(詳細についてはブログ記事を参照してください)。コアを4つ持っていることを活かして、旧型では行えなかった共有メモリ型の並列プログラミング学習も行えます。

もし今からRaspberry Piスパコンを作るのであれば、 学習の目的に応じて、新型・旧型を使い分ければ良いと思います。

おわりに

本連載を読んでいただき、ありがとうございました。

第0回でも述べたように、 スパコンは私たちの生活をより豊かなものとするため必要不可欠なものです。 スパコンを自作するという本連載を通して、 スパコンの必要性や、高度な技術により支えられていることを感じていただき、 興味を持っていただけたなら幸いです。

最後に、スパコンについてより詳しく知りたい方のため、 本連載を書く際に参考にした書籍をご紹介します。

・絵でわかるスーパーコンピュータ 姫野龍太郎 著 ISBN: 978-4061547643
・スパコンプログラミング入門 片桐孝洋 著 ISBN: 978-4130624534
・スーパーコンピューターを20万円で創る 伊藤智義 著 ISBN: 978-4087203950
・ARMで学ぶ アセンブリ言語入門 出村成和 著 ISBN: 978-4863541269
・ルーター自作でわかるパケットの流れ 小俣光之 著 ISBN: 978-4774147451

 
 

著者プロフィール

西永俊文 1992年生まれの情報系大学生。一人前のエンジニアを目指して日々勉強中。PDA全盛期の影響から、組み込み機器が好き。今欲しい物は、現代の技術で作られた物理キーボード付きで4インチ以下の小さなスマートフォン。著書に『BareMetalで遊ぶ Raspberry Pi』がある。
〈四コママンガ:

OLYMPUS DIGITAL CAMERA
著者近影(手にしているのは今回製作したRasPiクラスター初号機)

連載目次

2014-12-01 RaspberryPiクラスタ製作記 第0回「スパコンって本当にスゴイの?」
2014-12-01 PHASE/0のインストール方法
2014-12-16 RaspberryPiクラスタ製作記 第1回「スパコンを作ろう」
2014-12-16 Raspberry Piの初期設定
2014-12-29 RaspberryPiクラスタ製作記 第2回「スパコンで遊ぼう」
2014-12-29 HPLのインストール方法
2015-02-11 RaspberryPiクラスタ製作記 第3回「並列プログラミング」
2015-02-11 MPIを用いた並列プログラミングの概要
2015-02-11 LU分解アルゴリズムのおさらい
2015-03-01 RaspberryPiクラスタ製作記 第4回「BareMetal並列計算」

RaspberryPiクラスタ製作記 第3回「並列プログラミング」

第3回 並列計算プログラムをつくろう

mang3

前回までは既存のプログラムを使い、並列化による性能上昇について調べました。これにより大雑把に「並列化するとプログラムの処理を高速化できる」ことはわかりましたが、具体的にどのようにすれば自分のプログラムを並列化して高速化できるのかは、まだわかりません。

そこで今回からは、連立一次方程式を解くプログラムの制作を通して、並列計算プログラムを1から自作する方法や、高速化のための手法について学んでいきたいと思います。

なぜ連立一次方程式を扱うかというと、それは理解が容易だからです。スパコンでは様々な計算が行われていますが、大抵これらの計算は非常に難しく、短時間で理解し、プログラムを作るのは大変です。しかし、連立一次方程式は小さいサイズなら筆算で簡単に解けるため、これを解くプログラムの理解も容易です。加えて、連立一次方程式をプログラムで解く方法は、数値計算の教科書や並列計算の入門本など、非常に多くの書籍で扱われています。 そのため、内容でわからないことが出た場合や、より詳しく知りたい場合などは、豊富な参考文献をあたって調べられます。

本記事も、片桐孝洋さんの『スパコンプログラミング入門 並列処理とMPIの学習』という本と、同じく片桐さんの講義資料を参考にしています。わからないことが出た際はこちらを参照してください。

LU分解法(1)
LU分解法(2)
LU分解法(3)
東京大学情報基盤センター 准教授 片桐孝洋

連立一次方程式をプログラムで解く

今回は「LU分解」を用いて、以下の簡単な連立一次方程式を解くプログラムを作っていきます

eq01
式01

一般に、問題を解くプログラムのことをソルバと呼び、今回の連立一次方程式を解くプログラムのことは連立一次方程式ソルバと呼びます。

LU分解のアルゴリズムはいくつかあります。 今回は片桐さんの資料でも並列化に向くと紹介されている、外積形式ガウス法を用いて計算を行います。

※外積形式ガウス法を用いたLU分解のアルゴリズムについては下記のページを参照してください
アルゴリズムのおさらい

行列の定義

係数行列A、変数ベクトルx、定数項b、ベクトルcはプログラムで次のように定義します。

// SIZE: 係数行列Aのサイズ
#define SIZE 3

double A[SIZE][SIZE] = {
    {4, 1, 1},
    {1, 3, 1},
    {2, 1, 5}
};

double b[SIZE]={9,10,19};

double c[SIZE];
double x[SIZE];

LU分解のコードを作る

最初に係数行列Aを下三角行列Lと上三角行列Uに分解するコードを作ります。

下三角行列Lと上三角行列Uは一般的に1つの行列で表されます。

eq02
式02

行列LUの要素は、LU分解より係数行列Aを更新しつつ作られるので、 係数行列Aに上書きして入れていきます。 そして、LU分解が完了すると、行列Aの中身はすべて行列LUの要素におきかわります。 その後の前進代入と後進代入は、行列Aを置き換えて作られた行列LUを用いて行います。

rasp4f1

外積形式ガウス法を用いたLU分解は、 行方向のループk(k=0,1,…,n=SIZE)と、 そこに入れ子になった以下の2つのループで行います。

・行列L(枢軸ベクトル)を作っていく列方向のループ……(1)
・行列Uを作り、残りの値を更新する行と列の2重ループ……(2)

です。

raspi4f2

この時、1つ目のループで計算されるベクトルのことを枢軸(すうじく)ベクトルと呼びます。
以上の処理をプログラムで表すと、次のようになります。

int n = SIZE;

// LU分解開始
for (k=0; k<n; k++) {
    // kで回す行方向ループ

    // (1)
    // 行列Lの処理
    double div_tmp = 1.0 / A[k][k];
    for (i=k+1; i<n; i++) {
        // iで回す行方向ループ

        // l_ik = a_ik / a_kk
        A[i][k] = A[i][k] * div_tmp;
    }

    // (2)
    // 行列Uの処理
    for (j=k+1; j<n; j++) {
        // jで回す行方向のループ

        // Tips: ここで値を変数に入れたほうが
        // ループ内で参照するよりメモリアクセス回数が減ってお得
        double a_jk = dtemp = A[j][k];

        for (i=k+1; i<n; i++) {
            // iで回す列方向のループ
            A[j][i] = A[j][i] - A[k][i] * a_jk;
        }
    }
}
// LU分解終了

前進代入

LU分解の後は下三角行列Lと定数項bの値より、Lc=bのベクトルcの値を求めます。

// 前進代入
// Lc = b
for(k=0; k<n; k++){
    // 行方向のループ
    c[k] = b[k];
    for(i=0; i<k; i++){
        // 列方向のループ
        c[k] -= A[k][i] * c[i];
    }
}

後進代入

前進代入で求めたベクトルcより、Ux = cを解いて変数ベクトルxを求めます。

// 後進代入
// Ux = c
// k = n - 1
x[n-1] = c[n-1] / A[n-1][n-1];
// k = (n - 2) ~ (0)
for(k = n - 2; k >= 0; k--){
    x[k] = c[k];
    for(i=k+1; i<n; i++){
        // 列方向のループ
        x[k] -= A[k][i] * x[i];
    }
    x[k] = x[k] / A[k][k];
}

完成

以上でLU分解を用いた連立一次方程式ソルバは完成です。 以下にコードを示します。

LU_normal.c

このコードを実行すると、式eq01の係数行列AをLU分解した結果と、変数ベクトルxの解が出力されます。

$ gcc -o LU_normal LU_normal.c

$ ./LU_normal
LU:
 4.00  1.00  1.00
 0.25  2.75  0.75
 0.50  0.18  4.36

x:
 1.00  2.00  3.00

LU分解コードを並列化しよう

ここからはMPIを用いて、連立一次方程式を並列に解くプログラムを作っていきます。 ただし全ての並列化を行うのは大変なので、計算に一番時間がかかるLU分解の部分のみを並列化します。

※MPIについてご存じない方は、以下のページを参照してください。
MPIを用いた並列プログラミングの概要

最初に完成したコードを示し、このコードを扱いながら説明を行います。 以下からダウンロードして下さい。

LU_mpi_3x3

担当範囲の計算

外積形式ガウス法を用いたLU分解の並列化を、今回は列方向に分割して行います。

raspi4f3

各プロセスはそれぞれ、以下のプログラムで計算した、 開始(start)から終了(end)までの列を担当し、計算します。 担当列の長さは(partition_size)で表します。

// 1つのプロセスが担当する列の長さ
int partition_size = n / numprocs;

// プロセスが担当する列の開始点
int start = myid * partition_size;

// プロセスが担当する列の終了点
int end = (myid + 1) * partition_size;

行列Lの計算と放送

raspi4f5

下三角行列L(枢軸ベクトル)の計算は、その列を担当するプロセスが行います。

枢軸ベクトルは、その他のプロセスが係数行列Aの更新に必要とするので、別のプロセスに放送します。 そのため、枢軸ベクトルを担当したプロセスは、 別途放送用の配列Lを用意し、更新した値を配列Lに入れて放送します。

この処理を行うのが、次のコードです。

int worker = k / partition_size;

// 行列Lの処理
if(worker == myid){
    // kが自分の担当範囲内の場合はLを処理
    double div_tmp = 1.0 / A[k][k];
    for (i=k+1; i<n; i++) {
        // iで回す行方向ループ

        // l_ik = a_ik / a_kk
        A[i][k] = A[i][k] * div_tmp;
        L[i] = A[i][k];
    }
}

// Lを放送
MPI_Bcast(&L[k], n - k, MPI_DOUBLE, worker, MPI_COMM_WORLD);

行列Uを計算する

枢軸ベクトルを受け取った各プロセスは、行列Uの計算と、行列Aの値の更新を行います。

raspi4f6

最初に行うのは、計算に使用する枢軸ベクトルの要素a_jkの取り出しです。

double a_jk = L[j];

次に自身の担当する列のみ計算するため、i-ループを次のように書き換えます。

if(k < start){
    i = start;
}
else{
    i = k+1;
}

for (; i < end; i++) {
// 行列Uの計算と更新処理

kがプロセスの担当列番号より小さい場合は、startからendまでループを回して計算します。

その他の場合はkが自分の担当範囲か、自分の担当範囲が既に終了しているので、k+1からendまでループを回して計算します。

以上でLU分解の処理は終了です。

結果を集計する

今回の並列化はLU分解のみとしているので、 LU分解後の処理はすべてプロセス0が行います。 そのため、LU分解後は各プロセスの計算した配列LUのデータを、プロセス0に集めます。

これを行うのが次のコードです。

/ 結果を集める
for(i = 0; i < n; i++){
    // 行方向のループ
    for(j = partition_size; j < n; j+=partition_size){
        // 列方向のループ

        int worker = j / partition_size;

        if(myid == 0){
            // process0 はその他のプロセスの担当部分を集める
            MPI_Recv(&A[i][j], partition_size, MPI_DOUBLE, worker, 0, MPI_COMM_WORLD, &istatus);
        }
        if(myid == worker){
            // process0 以外は、自分が担当した部分を送る
            MPI_Send(&A[i][j], partition_size, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
        }
        // すべてのプロセスがここに来るのを待つ
        MPI_Barrier(MPI_COMM_WORLD);
    }
}

性能測定

プログラムの並列化は少々大変でしたが、 このプログラムを複数ノードで並列に動かすことで、 高速に計算できるはずです。

係数行列のサイズは3程度では変化がわからないので、 サイズを192に拡大して実行しました。 その結果をグラフに示します。

graph01

PHASE/0のときのようにノード数が増えるごとに計算時間が減って……いませんね。 今回のプログラムはPHASE/0の時と真逆で、ノード数を増やすごとに処理時間が増えています。

どうやらLU分解を並列化しただけでは、並列化のコスト(主に通信のコスト)のほうが大きいため、逆に時間が増えてしまうようです。

残りの前進代入と後進代入も並列化すれば、ノード数に応じて処理時間が減っていく、期待通りの結果が得られるかもしれません。余力のある方は、最初に挙げた片桐さんの資料を参考に、並列化を行ってみてください。私は別の方法で高速化を図ってみます。

BareMetalで作る連立一次方程式ソルバ専用機

現代のOSは様々な機能を私たちに提供してくれています。マルチタスクシステムもその1つです。このシステムは同じマシンを複数の人が使うことを想定し、プロセスを人間が分からない程度の短時間で切り替え、あたかも複数のプログラムが同時に動いているかのように実行するシステムです。これのおかげで私たちは先程のLU分解のプログラムを動かしている間も、別マシンからSSHアクセスして作業することが可能になっています。

しかし、今回のように私1人が1つのプログラムだけを処理したい場合、マルチタスクシステムを含め、OSの便利な機能は必要のない、無駄な機能です。このマルチタスク処理を含め、OSの実行にかかる様々なコストを削減して高速化できないでしょうか。

私はこの目的を達成する方法に、心当たりがあります。著者紹介欄でもご紹介頂いているように、私は「Baremetalで遊ぶ Raspberry Pi」という電子書籍を書きました。

通常Raspberry Piは、電源投入後にSDカード上のプログラムを読み込み、 LinuxなどのOS(カーネル)を起動します。 ここで、SDカード上の設定を少し書き換えると、本来カーネルが起動するところを、代わりに自作のコードを実行できます。そうすればOSの実行にかかる様々なコストを削減して、計算を高速に行えるはずです。

raspi4f7

この自作のコードの作り方と、OSを用いない素のハードウェア(BareMetal)上で、 Raspberry Piを制御する方法を説明した電子書籍が「Baremetalで遊ぶ Raspberry Pi」です。

ここからは前半に作った連立一次方程式ソルバをBareMetal上に再実装して、高速化を測ります。これはいわば、Raspberry Piを1次方程式を解くための専用機として運用するようなものです。このような贅沢な使い方ができるのも、自作スパコンの特権です。

早速、製作を行っていきましょう。

クロスコンパイル開発環境を整える

プログラムのバイナリを動かす環境と、コンパイルを行う環境が異なるコンパイルのことを、クロスコンパイルといいます。 BareMetal開発を行うには専用のクロスコンパイル開発環境が必要なので、 これから用意していきます。

以降の作業をRaspberry Piで行うと、最低でもまる1日以上かかってしまいます。 なので性能の良い別のマシン上(以降、母艦)にクロスコンパイル環境を用意します。

環境を合わせるため、母艦のOSはDebian Linux(wheezy) x86_64版を用います。 他の環境を使われている方は、仮想マシン上にこの環境を用意して行ってください。

Crosstool-NGのインストール

クロスコンパイル環境作成は、Crosstool-NGというツールを使って行います。 以下のコマンドを実行してインストールしてください。

$ wget http://crosstool-ng.org/download/crosstool-ng/crosstool-ng-1.20.0.tar.bz2
$ tar jxvf crosstool-ng-1.20.0.tar.bz2
$ cd crosstool-ng-1.20.0
$ ./configure
$ make
$ sudo make install

クロスコンパイラのビルド

ここからはCrosstool-NGを使ってクロスコンパイル環境を作成します。 作成に必要な設定ファイルを以下からダウンロードして、母艦のホームディレクトリに保存してください。

hard_dot_config

ダウンロード完了後、 以下のコマンドを実行すると、 クロスコンパイル環境のビルドが始まります。 環境にもよりますが、ビルド完了までは大体30分以上かかります。 気長にお待ちいただければと思います。

$ mkdir -p ~/ct-ng-rpi/src
$ cd ~/ct-ng-rpi
$ ct-ng arm-unknown-eabi
$ mv ~/hard_dot_config .config
$ ct-ng build

パスの追加

作成の終わったクロスコンパイル環境は ~/cross/arm-hardfloat-eabi 以下に保存されます。 このツールを使えるよう、システムにパスを通します。 以下のコマンドを実行してください。

$ export PATH=$PATH:$HOME/cross/arm-hardfloat-eabi/bin
$ echo 'export PATH=$PATH:$HOME/cross/arm-hardfloat-eabi/bin' >> ~/.bashrc

これで開発環境は整いました。

プログラムのクロスコンパイルとインストール

BareMetal環境上ではOSのもつ資源が一切ないので、 初期状態ではprintfどころかC言語自体がまともに動きません。

この問題を解決するには自分で初期化コードを書く必要があります。しかし、この作業は今回の本質でないので省略して、著書のコードをベースに開発を行います (この内容について興味のある方は、著書を読んでいただけると嬉しいです)。

以下に著書のコードをベースとして、 今回作成したLU分解の連立一次方程式ソルバをMyLUsolver関数として実装したコードを示します。次のリンクからダウンロードして下さい。

LU_baremetal_normal.zip

クロスコンパイル

先ほどダウンロードしたコードをクロスコンパイルして、 BareMetal Raspberry Piで動くコードを生成します。

コードにはMakefileを用意してあるので、makeコマンドでコンパイルは完了します。

$ make

実際に動かすバイナリの名前は rpi-micon.imgです。

インストール

Raspberry Piは起動時に、第1パーティションの設定ファイル config.txt を読み込んで、各種初期設定を行います。 この設定ファイルを書き換えると、起動時に読み込むプログラムを書き換えられます。

ここからはその設定ファイルを書き換え、自作コードが動くように設定を行います。 Raspberry Piに刺していたSDカードを抜き、パソコンに接続してください。

img06

するとこのようにBOOT(第1パーティション)とRECOBERYの2つのパーティションが見えるはずです。

img07

BOOTパーティションを開くと、中にconfig.txtが見つかります。 このファイルの末尾に、以下の1行を追記してください。

kernel=rpi-micon.img

これでRaspberry Piは、第1パーティションのrpi-micon.imgを起動するようになります。

最後にrpi-micon.imgを母艦からSDカードに移せばインストールは完了です。

シリアル通信

BareMetal環境上の標準出力はRaspberry PiのUART(シリアル通信)に出ていて、信号レベルは3.3Vです。 これを母艦で読み取るため、USB-シリアル変換ケーブルとRaspberry Piを接続します。

※シリアルケーブルはFTDI社のケーブル(秋月電子)がお勧めです。

Raspberry PiのUARTはP1コネクタのGPIO14番からTX、GPIO1510番からRXが出ています。

img08

上記の図と、USB-シリアル変換ケーブルの仕様を参照して
・Raspberry PiのTXをUSB-シリアル変換ケーブルのRX
・Raspberry PiのRXをUSB-シリアル変換ケーブルのTX
・Raspberry PiのGNDをUSB-シリアル変換ケーブルのGND
に接続してください。

シリアル通信を行うプログラムは、Windowsの方はTeraTerm、Macの方はCoolTermをおすすめします。ボーレートは115200bpsに設定してください。

BareMetal版の性能測定

SDカードをRaspberry Piに刺し電源を起動するとプログラムが動作し、連立一次方程式の解と実行時間を出力します。

最初に比較のため、Linux上で動かした連立一次方程式ソルバの結果を示します。 コードは以下からダウンロードして下さい。

LU_normal_192x192.zip

$ gcc -O2 -o LU_normal_192x192 LU_normal_192x192.c
$ ./LU_normal_192x192

exec time: 0.380212(sec)

結果は380ミリ秒でした。 これに比べてBareMetalでは、はたしてどれだけ性能が向上するでしょうか。 以下に結果を示します。

exec time: 458917(us)

459ミリ秒とLinux版に対して20%ほど余計に時間がかかっています。

……忘れていました。 OSが無いということは、ただ開発が不便なだけではないのでした。起動直後のCPUは、ハードウェア浮動小数点演算器(VFP)やキャッシュなど、様々な機能がOFFになっています。 そのため起動時にはこれらの機能を適切に設定することで、ハードウェアの性能を存分に発揮できる状態にしてあげなければいけません。 それを行うコードもOSのカーネルに組み込まれています。

私は今回、VFPの設定は行いましたが、キャッシュの設定を行っていませんでした。そのためどうやら、キャッシュメモリの設定が行われているLinuxの方は、命令かデータがキャッシュに残っており、その差が実行時間の差に現れたのだと考えられます。

リベンジ——VFPを効率的に用いる最適化

このままでは悔しいので、BareMetal版コードの方をいじって高速化してみます。

Raspberry PiのCPUにはVFP(v2)というハードウェア浮動小数点演算器がついており、 これを使うと浮動小数点演算を高速に処理できます。早速やってみましょう。

通常、ARMのCPUで計算を行うには、CPUの持つ汎用レジスタ(r0,r1,...,r12)を用います。 しかし、VFPはCPUとは独立した装置なので、VFPの計算はVFPが持つレジスタ間のみで行われます。 そのため計算を行うには、値をメモリからVFPのレジスタにコピーし、結果をメモリにストアする作業が行われます。

一般的に、メモリのロード・ストアや、異なる装置間の値交換は、 レジスタのみを用いた演算に比べ低速です。 しかし、この遅い処理を一度にまとめて行うと、パイプラインが効いて速く処理できる場合があります。 この恩恵を最大限に受けられるようにしましょう。

ループのアンローリング

Raspberry PiのVFPは16本の64bitレジスタ(d0~d15)を持っているので、 最大16個の値を一度にまとめて処理して効率化できます。 これがどこに適応できるかというと、LU分解の行列Lを処理した以下の部分と、

// 行列Lの処理
double div_tmp = 1.0 / A[k][k];

/////////////////////////////////
////// 適応できる場所ここから //////
/////////////////////////////////
for (i=k+1; i<n; i++) {
    // iで回す行方向ループ

    // l_ik = a_ik / a_kk
    A[i][k] = A[i][k] * div_tmp;
}
/////////////////////////////////
////// 適応できる場所ここまで //////
/////////////////////////////////

並列化の時にも触った、行列Uを処理する以下の部分です。

// 行列Uの処理
for (j=k+1; j<n; j++) {
    // jで回す行方向のループ

    // Tips: ここで値を変数に入れたほうが
    // ループ内で参照するよりメモリアクセス回数が減ってお得
    double a_jk = A[j][k];

    /////////////////////////////////
    ////// 適応できる場所ここから //////
    /////////////////////////////////
    for (i=k+1; i<n; i++) {
        // iで回す列方向のループ
        A[j][i] = A[j][i] - A[k][i] * a_jk;
    }
    /////////////////////////////////
    ////// 適応できる場所ここまで //////
    /////////////////////////////////
}

ここを一度にまとめて行うよう、次のように書き換えます。

/ 行列Lの処理
double div_tmp = 1.0 / A[k][k];
for (i=k+1; i<n-12; i+=12) {
    // iで回す行方向ループ

    double d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15;

    // l_ik = a_ik / a_kk
    d4 = A[i][k];
    d5 = A[i+1][k];
    d6 = A[i+2][k];
    d7 = A[i+3][k];
    d8 = A[i+4][k];
    d9 = A[i+5][k];
    d10 = A[i+6][k];
    d11 = A[i+7][k];
    d12 = A[i+8][k];
    d13 = A[i+9][k];
    d14 = A[i+10][k];
    d15 = A[i+11][k];

    d4 *= div_tmp;
    d5 *= div_tmp;
    d6 *= div_tmp;
    d7 *= div_tmp;
    d8 *= div_tmp;
    d9 *= div_tmp;
    d10 *= div_tmp;
    d11 *= div_tmp;
    d12 *= div_tmp;
    d13 *= div_tmp;
    d14 *= div_tmp;
    d15 *= div_tmp;

    A[i][k] = d4;
    A[i+1][k] = d5;
    A[i+2][k] = d6;
    A[i+3][k] = d7;
    A[i+4][k] = d8;
    A[i+5][k] = d9;
    A[i+6][k] = d10;
    A[i+7][k] = d11;
    A[i+8][k] = d12;
    A[i+9][k] = d13;
    A[i+10][k] = d14;
    A[i+11][k] = d15;
}
for(;i<n;i++){
    A[i][k] = A[i][k] * div_tmp;
}
// 行列Uの処理
for (j=k+1; j<n; j++) {
    // jで回す行方向のループ

    // Tips: ここで値を変数に入れたほうが
    // ループ内で参照するよりメモリアクセス回数が減ってお得
    double a_jk = A[j][k];

    for (i=k+1; i<n-6; i+=6) {
        // iで回す列方向ループ
        double d4,d5,d6,d7,d8,d9;

        // l_ik = a_ik / a_kk
        d4 = A[k][i+0];
        d5 = A[k][i+1];
        d6 = A[k][i+2];
        d7 = A[k][i+3];
        d8 = A[k][i+4];
        d9 = A[k][i+5];

        d4 *= a_jk;
        d5 *= a_jk;
        d6 *= a_jk;
        d7 *= a_jk;
        d8 *= a_jk;
        d9 *= a_jk;

        A[j][i+0] -= d4;
        A[j][i+1] -= d5;
        A[j][i+2] -= d6;
        A[j][i+3] -= d7;
        A[j][i+4] -= d8;
        A[j][i+5] -= d9;
    }
    for(;i<n;i++){
        A[j][i] = A[j][i] - A[k][i] * a_jk;
    }
}

このようにループを展開して高速化を行う手法を、ループのアンローリングといいます。

行列Lで行っていた計算は、

A[i][k] = A[i][k] * div_tmp;  

このようなdiv_tmpと行列Aの要素1つの単純な掛け算なので、使用するVFPのレジスタ16以下に収まるよう、12回分(1+12 < 16)を一度に行うよう展開しています。 行列Uの計算と更新は、

A[j][i] = A[j][i] - A[k][i] * a_jk;

このようにa_jkと行列Aの要素を掛けた後、行列Aの別の要素へ引き算を行っています。 これはa_jkに1つ、行列Aの要素に2つレジスタを使うので、6回分(1+2*6 < 16)を一気に行うよう展開しています。

ループのアンローリング効果の確認

ループのアンローリングによる効果は、逆アセンブル結果を見れば明らかです。 以下にアンローリング後のコードを示すので、ダウンロードしてmakeを行ってください。

LU_baremetal_unrolling.zip

逆アセンブル結果は rpi-micon.disas に出力されます。このファイルを、アンローリング前と後で見比べてみてみます。

行列Lを処理していた部分の、アンローリング前のコードを示します。

809c:   ed937b00    vldr    d7, [r3]
80a0:   ee277b06    vmul.f64    d7, d7, d6
80a4:   ed837b00    vstr    d7, [r3]
80a8:   e2833c06    add r3, r3, #1536   ; 0x600
80ac:   e15c0003    cmp ip, r3
80b0:   1afffff9    bne 809c <myLUsolver+0x34>

このコードはループを使って処理しているので、1変数ずつvldrでロード、vmul.f64で乗算、vstrでストアを繰り返しています。 しかしアンローリングを行うと、次のようにロードとストアが一列に並ぶコードになります。

8124:   ed92ab00    vldr    d10, [r2]
8128:   ed919b00    vldr    d9, [r1]
812c:   ed9b8b00    vldr    d8, [fp]
8130:   ed9a0b00    vldr    d0, [sl]
8134:   ed991b00    vldr    d1, [r9]
8138:   ed982b00    vldr    d2, [r8]
813c:   ed953b00    vldr    d3, [r5]
8140:   ed944b00    vldr    d4, [r4]
8144:   ed9c5b00    vldr    d5, [ip]
8148:   ed93bb00    vldr    d11, [r3]
814c:   ed97cb00    vldr    d12, [r7]
8150:   ed906b00    vldr    d6, [r0]
8154:   ee2aab07    vmul.f64    d10, d10, d7
8158:   ee299b07    vmul.f64    d9, d9, d7
815c:   ee288b07    vmul.f64    d8, d8, d7
8160:   ee200b07    vmul.f64    d0, d0, d7
8164:   ee211b07    vmul.f64    d1, d1, d7
8168:   ee222b07    vmul.f64    d2, d2, d7
816c:   ee233b07    vmul.f64    d3, d3, d7
8170:   ee244b07    vmul.f64    d4, d4, d7
8174:   ee255b07    vmul.f64    d5, d5, d7
8178:   ee2bbb07    vmul.f64    d11, d11, d7
817c:   ee2ccb07    vmul.f64    d12, d12, d7
8180:   ee266b07    vmul.f64    d6, d6, d7
8184:   e59d3000    ldr r3, [sp]
8188:   e283300c    add r3, r3, #12
818c:   e58d3000    str r3, [sp]
8190:   e35300b3    cmp r3, #179    ; 0xb3
8194:   e59d3008    ldr r3, [sp, #8]
8198:   ed82ab00    vstr    d10, [r2]
819c:   ed819b00    vstr    d9, [r1]
81a0:   ed8b8b00    vstr    d8, [fp]
81a4:   ed8a0b00    vstr    d0, [sl]
81a8:   ed891b00    vstr    d1, [r9]
81ac:   ed882b00    vstr    d2, [r8]
81b0:   ed853b00    vstr    d3, [r5]
81b4:   ed844b00    vstr    d4, [r4]
81b8:   ed8c5b00    vstr    d5, [ip]
81bc:   ed83bb00    vstr    d11, [r3]
81c0:   e2833b12    add r3, r3, #18432  ; 0x4800
81c4:   ed87cb00    vstr    d12, [r7]
81c8:   e2822b12    add r2, r2, #18432  ; 0x4800
81cc:   ed806b00    vstr    d6, [r0]
81d0:   e2811b12    add r1, r1, #18432  ; 0x4800
81d4:   e28bbb12    add fp, fp, #18432  ; 0x4800
81d8:   e28aab12    add sl, sl, #18432  ; 0x4800
81dc:   e2899b12    add r9, r9, #18432  ; 0x4800
81e0:   e2888b12    add r8, r8, #18432  ; 0x4800
81e4:   e2855b12    add r5, r5, #18432  ; 0x4800
81e8:   e2844b12    add r4, r4, #18432  ; 0x4800
81ec:   e28ccb12    add ip, ip, #18432  ; 0x4800
81f0:   e58d3008    str r3, [sp, #8]
81f4:   e2877b12    add r7, r7, #18432  ; 0x4800
81f8:   e2800b12    add r0, r0, #18432  ; 0x4800
81fc:   daffffc8    ble 8124 <myLUsolver+0xbc>

美しいですね。

連続してロードストアや演算を行っているので、パイプラインが効きやすいコードになっているはずです。

最適化結果の確認

ではこのプログラムを実行して、どれだけ最適化できたか確認してみましょう。 以下に結果を示します。

exec time: 417205(us)

アンローリングすると、その分コードサイズは増えてしまいますが、 命令の実効効率が良くなります。 その証拠に、アンローリング後のプログラムの実行時間は約40msも短くなりました。

しかし、Linux上で動かした結果には後一歩及びませんでした。リベンジ失敗です。

BareMetal版が勝つためには……

BareMetal版がLinuxで動かした結果に勝つためには、 すくなくともLinuxで行っているハードウェアの設定を、BareMetalでも行う必要があるでしょう。 しかし、その設定を行うため、CPUの資料やコードを読むコストを考えると、BareMetalでやるより、Linuxカーネルの設定を修正したほうが速いでしょう。

BareMetal版が優っている点としては、電源投入後すぐに実行され、結果が出ることでしょうか。 なので、電源投入から解が出るまでの時間で競えば、BareMetal版の勝利です。

今回のまとめ

・並列プログラミングの勉強には連立一次方程式ソルバがおすすめ
・並列化しても高速にならないどころか、逆に遅くなることもある
・OSの資源は偉大
・ループのアンローリングを適切に行うと、処理を高速化できる

次回は通信コストを削減して、再度、並列処理に挑戦する予定です

著者プロフィール

西永俊文 1992年生まれの情報系大学生。一人前のエンジニアを目指して日々勉強中。PDA全盛期の影響から、組み込み機器が好き。今欲しい物は、現代の技術で作られた物理キーボード付きで4インチ以下の小さなスマートフォン。著書に『BareMetalで遊ぶ Raspberry Pi』がある。
〈四コママンガ:

OLYMPUS DIGITAL CAMERA
著者近影(手にしているのは今回製作したRasPiクラスター初号機)

連載目次

2014-12-01 RaspberryPiクラスタ製作記 第0回「スパコンって本当にスゴイの?」
2014-12-01 PHASE/0のインストール方法
2014-12-16 RaspberryPiクラスタ製作記 第1回「スパコンを作ろう」
2014-12-16 Raspberry Piの初期設定
2014-12-29 RaspberryPiクラスタ製作記 第2回「スパコンで遊ぼう」
2014-12-29 HPLのインストール方法
2015-02-11 RaspberryPiクラスタ製作記 第3回「並列プログラミング」
2015-02-11 MPIを用いた並列プログラミングの概要
2015-02-11 LU分解アルゴリズムのおさらい
2015-03-01 RaspberryPiクラスタ製作記 第4回「BareMetal並列計算」


↑前のページ次のページ↓