計算工学ナビ

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

サイト内検索

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並列計算」

RaspberryPiクラスタ製作記 第2回「スパコンで遊ぼう」

第2回 Raspberry Piスパコンで遊ぼう

自作したばかりのスーパコンピュータ上で、 第0回で動かした第一原理電子状態計算ソフトウェアPHASE/0を動かして遊んでみましょう。

この遊びの目的は以下の2つです。

・自作スパコンでも「京」で動いている計算プログラムが動くことを確認する
・並列化により処理がどれだけ高速化されるかを確認する

「並列化により処理がどれだけ高速化されるか」に対して、まずあえてもっとも単純に予想すると「並列計算ノード数を増やすと、それに比例して性能も上昇する」となります。ノード数を倍にするということは計算に使えるリソースが倍になるということなので、ノード数に比例して処理速度も上がっていくのではないでしょうか。この予想が正しいかどうかを、これから検証していきます。

PHASE/0のインストール

今回使用するRaspberry PiはARMアーキテクチャのCPUを搭載しているため、第0回でx86マシン用にコンパイルしたバイナリを実行することができません。 そこで再度ソースコードからのコンパイルを行います。

最初にPAHSE/0のコンパイルと動作に必要なパッケージをインストールします。 以下のコマンドを実行してください。

$ sudo aptitude -y install gfortran mpich2 libc6-dev

PHASE/0のコンパイルは基本的に ./install.sh を実行するだけですが、こちらで試したところ、今回使用したバージョンのインストールスクリプトはRaspbianでは正常に動作しませんでした。そこで、このインストールスクリプトが正常動作する環境(MacOSX)で動作させてコンパイルに必要なMakefileを生成し、これをRaspberry Pi上に持ち込むことでコンパイルを可能としました。このMakefileを以下のURLより取得してください。

Makefile

取得したMakefileを phase0_2014.03/src_phase/ 以下に置いて、そのディレクトリ上で make コマンドを実行するとコンパイルが開始されます。エラーなくコンパイルが終了すれば、PHASE/0のインストールは完了です。

計算打ち切り時間の書き換え

PHASE/0の入力ファイル中には、計算をどれだけで打ち切るか(cpumax)の設定が書かれています。今回使用するRaspberry Piは第0回で使用したマシンに比べて非常に性能が低いため、デフォルトの打ち切り時間では計算が終了しません。そのためこの設定を修正して、打ち切り時間を延長します。

samples/Si8/input_scf_Si8.data を開いて、以下の3600secとなっている部分を1dayに修正してください。

Control{
		# 修正前
        cpumax = 3600 sec     ! {sec|min|hour|day}
		# 修正後
        cpumax = 1 day     ! {sec|min|hour|day}
        condition=initial
}

MPIノードの設定

PHASE/0のプログラムはMPI方式でデータをやり取りする、mpichというアプリケーションを用いて並列計算を行います。 このmpichは、自分のマシンと一緒に計算を行うマシン(ノード)のホストネーム(またはIPアドレス)を設定ファイルに記述し、それを実行時に読み込ませることで並列計算が可能となります。 以下にその設定ファイルを記すので ~/mpdhosts に保存してください。

rpicluster01
rpicluster02
rpicluster03
rpicluster04
rpicluster05
rpicluster06
rpicluster07
rpicluster08
rpicluster09
rpicluster10
rpicluster11
rpicluster12
rpicluster13
rpicluster14
rpicluster15
rpicluster16

処理時間の計測

インストールの終わったPHASE/0を使い、 ノード数を1,4,9,16と変化させて、その処理時間を計測します。

1ノード1並列

最初は動作確認のため1並列で動かしてみましょう。 以下のコマンドで計算を開始してください。

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=1 ../../bin/phase -ne=1 -nk=1
real    32m59.518s
user    29m53.770s
sys     3m2.960s

実行時間は実時間で約33分と、第0回と比べて非常に時間のかかる結果となりました。 今回使用したRaspberry Piはクロックが700MHzの32bitマシンなので、妥当なスコアでしょう。

4ノード4並列

次にRaspberry Piクラスタマシンのノードを4台使用し、計算を4並列で回してみます。 予想のとおりであれば、CPUの数が4倍になるので実行時間も1/4になるはずですが、果たしてどうなるでしょうか。

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=4 ../../bin/phase -ne=2 -nk=2
real    16m14.437s
user    13m50.940s
sys     2m20.780s

実行時間は実時間で約16分、1並列の約半分となりました。 1/4とまでは行きませんでしたが、1/2でもまずまずのスコアでしょう。

この調子で9,16並列で計算を行い、その性能の上昇を見ていきます。

9ノード9並列

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=9 ../../bin/phase -ne=3 -nk=3
real    11m22.431s
user    9m41.020s
sys     1m39.270s

16ノード16並列

pi@rpicluster01 ~/phase0_2014.02_rpi/samples/Si8 $ time mpirun 
-f ~/mpihosts -np=16 ../../bin/phase -ne=4 -nk=4
real    9m25.834s
user    7m11.440s
sys     2m11.720s

計算結果の考察

1,4,9,16並列の計算時間の結果をまとめたグラフを示します。

exec_time

予想ではノード数を上げると計算に使えるCPUの数が増えるので、ノード数に比例して実行時間は減っていくはずでした。 しかし実際は、4並列の時の実行時間は一気に約1/2となったものの、 その後の9並列では1並列に比べ約1/3と伸びが悪くなり、16並列ではより一層伸びが悪くなっています。

一体どうして、計算時間はノード数に比例して下がっていかなかったのでしょうか。

アムダールの法則

計算時間がノード数に比例して下がらない理由ついて、第1回と同じように夏休みの宿題にたとえて説明します。 第1回では、夏休みの宿題として計算問題を100題だされたとしても、 100人の協力者が1問づつ解いて最後に答えを共有すれば、一瞬で解くことができると言っていました。 この計算問題を100人で解く手順を改めて示すと、以下の5手順に分けられます。

1.計算を100個に分割する
2.分けた計算を担当する協力者に解いてもらうよう指示を出す
3.それぞれ協力者が計算する
4.結果を集める
5.結果を問題集に書き込む

こうしてみると、100人の協力者が同時に作業を行える部分、つまり並列処理できる部分は手順3.だけです。 それ以外の部分は指示を出す人が1人で行う必要がある部分、つまり並列処理することができない部分です。

この例からわかるように、処理の中には並列に処理することにより高速化できる部分と、そうでない部分が存在します。 そのためノード数を上げることによる処理時間の減少は、一定のラインを超えると並列処理により高速化できない部分がネックとなり、処理時間の減少の伸びが非常に悪くなります。これが今回私たちの実験で起こった現象です。

このノード数(プロセッサ)と処理速度の関係には「アムダールの法則」という法則があります。 アムダールの法則は、並列化できない部分の割合をs、プロセッサ数をnとすると、 次の式で処理速度の向上比を計算することができます。

eq0

実際にこの式を使い、s=0.001、n=1~30000で計算した結果をグラフに示します。

exec_time

性能評価とコストの比較

コンピュータの性能には、理論性能と実行性能の2つがあります。

理論性能 ハードウェアの性能から単純に計算した最大性能。例えば700MHzで動作するRaspberry Piは、1クロックで1回浮動小数点演算が行えるので、理論性能700MFLOPSとなる。

実効性能 何らかのベンチマークソフトを動かして得た性能。

理論性能はあくまで理論的な最大性能なので、実効性能とは離れた結果となります。 先ほど私たちの自作スパコン上でPHASE0を動かした結果からもわかると思います。 しかし、具体的にどれだけの性能を持っているかは、まだ調べていませんでした。

性能の他にもうひとつ気になるのは、そのお値段です。現代はコストパフォーマンスを重視する傾向にあります。そのため、どれだけスパコンの性能が良くとも製作費やランニングコストが高すぎては実用的ではありません。

ところで、Raspberry PiのCPUアーキテクチャはARMというものでした。 このARMのCPUは性能に対して消費電力が低いという特徴があり、 身近なところでは携帯ゲーム機やスマートフォンで使用されています。 そのためもしかすると、消費電力――つまりランニングコストの面では、スパコン「京」に勝つことができるのではないでしょうか。

ここからはRaspberry Piスパコンの実効性能とそのコストを調べ、最後にスパコン「京」と比較してみましょう。

実効性能計測

Linpackは、コンピュータの実効性能の計測(ベンチマーク)を行うプログラムのひとつです。 TOP500ではこれを並列化した「HPL」(High-Performance Linpack)を用いて、ベンチマークスコアを測っています。 私たちのRaspberry Piクラスタマシンも、このHPLを用いて実効性能を計ってみましょう。

HPLのインストールについては、別のページで説明します。

HPLは処理する問題サイズ、ブロックサイズ等のパラメータ変更することで最適化を行い、 スコアを上げることができます。 しかし最適なパラメータを探すのは非常に大変なので、今回は私の試した中で一番よいスコアの得られた設定を利用します。

以下のファイルをダウンロードし、~/hpl/Linux_ARM_CBLAS/以下へ上書きコピーしてください。

HPL.dat

HPLの実行

実行方法はPHASE0とほとんどおなじですが、 こちらはノードの指定にホストネームを使うとエラーが発生するので、 IPアドレスを使ってノードを指定します。 本環境でのrpicluster01-16のIPアドレスは 192.168.240.101-116なので、16ノード計算を行うため以下のコマンドを実行します。

$ cd ~/hpl/bin/Linux_ARM_CBLAS
$ mpirun -np 16 -host 192.168.240.101,192.168.240.102, \
192.168.240.103,192.168.240.104, \
192.168.240.105,192.168.240.106, \
192.168.240.107,192.168.240.108, \
192.168.240.109,192.168.240.110, \
192.168.240.111,192.168.240.112, \
192.168.240.113,192.168.240.114, \
192.168.240.115,192.168.240.116  \
./xhpl

ベンチマーク結果

=================================================================
T/V                N    NB     P     Q               Time                 Gflops
-----------------------------------------------------------------
WR11C2R4       25600    96     4     4            4211.62              2.656e+00

HPLの結果は2.656GFLOPSとなりました。このスコアは、今から約30年前の1987年に発表されたスパコン「HITAC S-810」(最大3GFLOPS)に匹敵するスコアです。
 参考URL: http://museum.ipsj.or.jp/computer/super/index.html
当時は巨大なシステムで実現していた性能が、今ではこんな小さなボード16枚で得られてしまいます。 私は当時まだ生まれていませんが、技術の進歩が肌で感じられました。

ここからは自作Raspberry Piクラスタマシンの制作費とランニングコストを調べていきます。

製作費は11万円

今回製作に使用した部品とそのお値段の一覧を示します。 なお、製作の人件費はないものとします。

部品 単価(円)
Raspbebrry Pi TypeB+ 16 3,900
micro SDカード 8GB 16 800
USB-MicroBケーブル 16 600
2ポートUSB充電器 8 800
LANケーブル 16 350
24ポート スイッチングハブ 1 7,800
9ポート電源タップ 2 1,000
M2.6連結スペーサーネジ 72 10
小型メタルラック 1 600
結束バンド 1 100
ファンシー金網 1 100
金網用フック付き小物いれ 1 100

以上を合計すると、製作にかかったお金の合計は約11万円でした。

電気代(ランニングコスト)

USB充電器(Raspberry Piの外部電源)とUSBケーブルの間に簡易テスターを仕込んで、Raspberry Piクラスタマシンの消費電力を調べます。

ampere

結果は1台あたり200mA。Raspberry Piには負荷に応じてクロックを変えるような機能は無いため、ピーク時もこの値です。よって1台あたりの消費電力は 5(V) x 0.2(A) = 1(W)、16台でも16Wです。 その他、スイッチングハブやUSB充電器のロスを含めると、このクラスタマシンの総消費電力は大体20Wとなります。

この消費電力から電気料金を計算すると、1日あたり約10円となります。製作費同様、メンテナンスにかかる人件費は0円とすると、 ランニングコストは1日あたり10円、1年動かし続けても3560円です。

京との比較

比較のため、スパコン「京」のコストについても調べてみましょう。京のコストについての資料としては、平成25年1月に文部科学省からの資料があります。この資料によると、京の制作費は合計793億円(施設建設費用は含めず)、運用費は平成24年のデータで1年あたり97億円です。

ここまで調べたRaspberry Piクラスタと京の制作費を表に示します。

実効性能 制作費(円) 運用費(円/年)
10PFLOPS 793億 97億
RasPiクラスタ 2.656GFLOPS 11万 3560

これだけではわかりづらいので、 性能をコスト(製作費+運用費)で割った、1FLOPSあたりのコストを計算して比較します。

cost_graph

もしかするとコストの面では勝てるのではと考えていましたが、 仮に100年運用したとしても京を追い抜くことはできませんでした。技術の粋を集めて作られた京は、個人のレベルでは到底届かないものであることを再認識することが出来ました。

まとめ

今回はRaspberry Piスパコン制作のソフトウェアの整備と、 その上でPAHSE0の動作を行い、ノード数と計算時間の関係について調べてみました。 この内容を以下に簡単にまとめます。

・ノード数増加による処理速度向上は一定のところで伸び悩む
・アムダールの法則はプロセッサ数と性能向上比の関係を教えてくれる
・京のコストパフォーマンスは自作Raspberry Piスパコンと比べると圧倒的に良い

著者プロフィール

西永俊文 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並列計算」

12月5日(金)第5回分野4次世代ものづくりシンポジウム 開催報告

平成26年12月5日(金)、兵庫県神戸市にある(独)理化学研究所計算科学研究機構にて、文部科学省HPCI戦略プログラム分野4「次世代ものづくり」による、第5回「分野4 次世代ものづくり」シンポジウム−エクサスケールを見据えたものづくりシミュレーション−が開催されました。

当日は分野4各課題の研究開発成果進捗、産業界での利用事例等、計14件の発表がありました。多数の方々にご参加頂き、各セッション毎に活発な意見交換がなされました。

シンポジウムのプログラムについては以下のページを御覧ください。

http://www.ciss.iis.u-tokyo.ac.jp/supercomputer/event/event.php?id=55

DSCF4357

分野4統括責任者 東京大学生産技術研究所 加藤千幸教授による分野4全体進捗状況の説明

会場の様子

会場の様子


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