月別アーカイブ: 2013年4月

[C++11対応記念] clang/LLVMのインストールとGMPを使って円周率ベンチ。

Linuxの話です。

続編(?)はこちら: [フラクタル] マンデルブロー集合 10の200乗倍への挑戦

先日、注目のコンパイラclang/LLVMがC++11に完全対応したというニュースがあったので、せっかくの機会なのでどんな具合なのか試してみました。

円周率ベンチはだいぶ下です。結論を言ってしまうと、速度はそこまで大きく変わりませんでした。<

まずはclang/LLVMのインストール。更新も早いし、基本的には本家ページのガイドラインに従うのがよいと思います。C++11に対応したのはつい最近なので、ソースコードをリポジトリから入手します。

% mkdir llvm
% svn co http://llvm.org/svn/llvm-project/llvm/trunk llvm
% cd llvm/tools
% svn co http://llvm.org/svn/llvm-project/cfe/trunk clang
% cd ../projects
% svn co http://llvm.org/svn/llvm-project/compiler-rt/trunk compiler-rt

コンパイルにはコンパイラやgmakeなどいくつか基本的なソフトが必要になります。本家ページにも書いてありますが、だいたいのLinuxディストリビューションにはあらかた対応しているかと思います。
ただしGCC4.3以前やARMでの4.6以前だとコンパイルにいろいろ問題があるようです。

気分と環境が整ったらコンパイルとインストールをします。

% ./configure --enable-optimized --enable-jit --enable-targets=x86_64 --prefix=${HOME}/どこどこ
% gmake -j コア数
% make check
% make install

オプションは適当に。–prefixオプションはなくても構いませんが、今回とは別にclangを既にインストールしている、root権限がないのでホームディレクトリ以下にインストールしたい、などの場合には指定しておきましょう。
その場合はインストールディレクトリにパスを通しておくことをおすすめします。
gmake -j の後は、CPUのコア数に応じて数字を指定しておくとそれなりに早くなります。
手元の環境では10分もあれば終わりました。

GCCのインストールがあれこれ依存してて泣きたくなるのに比べたら、LLVMのインストールは非常にシンプルです。この辺りはさすが後発。

test1.cppにHello Worldを書いて確認。しかもあえて間違えて。

#include <iostream>

int main(int argc, char** argv){
  cout << "Hello, World!" << std::endl;
return 0;
}

そしてコンパイル。

% clang++ test1.cpp -o test1_clang
test1.cpp:4:2: error: use of undeclared identifier 'cout'; did you mean 'std::cout'?
        cout << "Hello, World!" << std::endl;
        ^~~~
        std::cout
/***/include/c++/4.7.2/iostream:62:18: note: 'std::cout' declared here
  extern ostream cout;          /// Linked to standard output
                 ^

なかなか賑やかでフレンドリーなエラーが出力されました。(といってもこの程度ならg++でも比較的穏やかなエラーを返してくれますが。)

coutをstd::coutに書きなおして実行。無事Hello Worldです。

インストールだけで満足してもつまらないので、円周率の計算でベンチマークをとってみます。円周率の計算ソフトはいろいろありますが、今回は多倍長計算でGMP (The GNU Multiple Precision Arithmetic Library)を導入してGMP謹製の円周率計算ベンチを回してみます。
(なんでBSDライクなライセンスのclang/LLVMを準備したのにGNUプロジェクトのソフトを持ってくるのかという感もありますが、まあGMPはLGPLですし許してください。)

なお、遊びでなく本気のベンチマークが見たい方はPhoronixのページなどを参考にして下さい。
ベンチマークに円周率を使ったのは気分です。

GMPのサイトに行って最新バージョンを持ってきます。現時点では5.1.1が最新。GCCとclangでコンパイルします。

% mkdir gmp
% cd gmp
% wget http://ftp.jaist.ac.jp/pub/GNU/gmp/gmp-5.1.1.tar.bz2
% tar xvfj gmp-5.1.1.tar.bz2
% cd gmp-5.1.1
% ./configure --prefix=${HOME}/どこどこ/gcc-4.7 --enable-cxx CC=gcc CXX=g++
% make
% make check
% make install
% ./configure --prefix=${HOME}/どこどこ/clang-3.3 --enable-cxx CC=clang CXX=clang++
% make
% make check
% make install

GMPの円周率計算ベンチのページからソースコードを持ってきて、適当にpi.cなどの名前で保存します。
そしてコンパイル。-Iと-Lのパスは適当に読み替えてください。

% gcc pi.c -I../gmp/gcc-4.7/include -L../gmp/gcc-4.7/lib -lgmp -lm -O3 -o pi-gcc
% clang pi.c -I../gmp/clang-3.3/include -L../gmp/clang-3.3/lib -lgmp -lm -O3 -o pi-clang3.3
/tmp/pi-fDvkj6.o: In function `bs':
pi.c:(.text+0xc90): undefined reference to `fac_mul2'
pi.c:(.text+0xd89): undefined reference to `fac_mul2'
pi.c:(.text+0xdf8): undefined reference to `fac_mul_bp'
pi.c:(.text+0xeb7): undefined reference to `fac_mul_bp'
pi.c:(.text+0xed6): undefined reference to `fac_mul_bp'
clang: error: linker command failed with exit code 1 (use -v to see invocation)

あれ?clangでリンカエラーが出るぞ?

ソースを見てみるとinline指定されている関数が見つからないことになっているもよう。プリプロセッサの問題かと思って-Eオプションつけてみても無事残っている。
…と調べてみたら本家に説明がありました。曰く、「C99の仕様ではinline指定をした関数定義はインライン展開されたときにのみ有効になるため、inline指定のない定義は別に必要になります。」とのこと。つまり、inlineと書いてあってもインライン展開されるかどうかはコンパイラが自由に決めることであり(clangでは「ゆるい提案」扱いになるとのこと)、インライン展開しないと決めたときに関数未定義になるらしい。

知らない。知らなすぎる。
clangで起きるエラーは実は手元のソースコードの仕様違反であった(がコンパイラの拡張やバグで潜在化していた)ということがよくあるようですが、今回もそうでした。コンパイラを変えてみるというのもそれだけで勉強になります。
まあ、今まで自分が知っていたと思っていた世界は実はC/C++でなくgcc/g++であった、というのが改めてよくわかります。

気を取り直して再コンパイル。inline 関数をstaticにするとかオプションに-std=gnu89をつけるとかいろいろありますが、単純にinlineを外すのが簡単です。どうせinlineつけてもclang最適化のほうが偉いので。

計測タイム。環境はXeon E5-2690 x2 (Sandy Bridge, 2.9/3.8GHz 2CPU/16コア/32スレッド)、メモリ16Gbyte、Red Hat Linux 64bit。 Xeonでコアが多いけど今回はマルチコア未対応なので関係なし。あえて言うならメモリ転送速度が4チャネルx2でやたら早い。

円周率100,000,000桁:
gcc : 111.126sec
clang : 111.309sec

あまり変わらない。桁をもっと増やしてみる。

円周率1,000,000,000桁:
gcc : 1664.053sec
clang : 1667.985sec

やっぱりほとんど変わらない。まあ、GMPが相当最適化されていてコンパイラの入る余地も少ないのでしょう。
個人的にはGMPをコンパイルできただけで私はだいぶ満足です。あと、円周率10億桁なんてとんでもないものが即座に計算できる時代になったことを再認識して喜んでいます。(clang関係ないけど)

それにしてもclangの進化はすごいですね。オープンソースでC/C++のコンパイラを作るなんて既にデファクトスタンダードのGCCもあるしで当分微妙なままなんじゃないかと思っていましたが、大間違いでした。セルフホスティングになってからあれよあれよという間にC++11完全対応まで来てしまいました。

clangはコンパイラとしてだけではなく、統合開発環境との連携も目玉のひとつです。今後Eclipseなど既存のIDEと深く関わっていけばコーディングがより気楽になるのでしょうか。

私も手元の遊びで書くC++はclangを使ってもいいかなと思っています。何よりエラーが見やすくて助かります。

[作ってみた] [プログラミング] 剛体3重振り子

3件の返信

以前ここで書きかけた剛体振り子のプログラムですが、いろいろと紛失してしまったので新しく作りました。

ソフトのダウンロードはこちらから pendulums

ソースコードはgithubからどうぞ https://github.com/So-Takamoto/triple_pendulums/
手元の環境はVisual Studio 2010です。

(以前の記事: [昔話] 振り子のシミュレーションをした話)

例によってC++, OpenGL+freeglutで作成しています。
キーの1,2,3で微分方程式を解く方法を変更でき、iで初期化、lで軌跡の表示切り替え、qで終了です。

質量や重力加速度、初期条件、各振り子の長さを変えたい場合はソースコードをいじってみてください。
64bit版でコンパイルを行いたい場合はfreeglutのライブラリを64bit版に変更してください。今回は軽いのであまり必要はないかと思いますが。

例によって動画を用意しておきます。

なんでこんなシミュレーションをしようと思ったかの話は以前書いた(早い話がこういうおもちゃが大学に置いてあった)ので、今回は原理とか内容の紹介メインにしたいと思います。

今回扱う系は平面上の棒3個なので、そのままで言えば3自由度×3個=9自由度の系なのですが、各ジョイントが2自由度分の拘束を持つため結局本来の自由度は3になります。

ひとつのやり方は、これらのジョイントにかかる力を陽に扱う方法です。利点は自由度を実座標と対応付けられるため意味がわかりやすいことと、比較的機械的な処理で立式できることです。
これらジョイントにかかる力は幾何学的な関係から導出することもできますし、今回のような単純なリンク機構であればエンドエフェクタ(一番先の振り子)から芋づる式に数値解析の範囲で導出することもできます。ロボットハンドなどの制御ではこの辺りの話が動力学計算としてたくさん出てきます。

一方で今回のような系ならそこまで面倒な話をしなくても、解析的に運動方程式を求めることができます。こういう系ではラグランジュ形式を出発点にすると便利です。ジョイントに働く力を陽に扱わずに、自由度の数だけ運動方程式を考えればいいからです。

ラグランジュの運動方程式の出発点は、ラグランジアンL(=T-U)の表現です。運動エネルギーTは重心の並進運動と回転運動から、位置エネルギーUは重心の位置から求まります。

3つの剛体を根本から順に0,1,2と置いて(0から始めるのはC言語使いの癖です)、それぞれ棒の角度を $\theta_{i}$ 、長さを $l_{i}$ 、質量を $m_{i}$ 、根元側のジョイントからの重心の相対位置を ${\bf x}_{ri}$ と置きます。3自由度はここでは3つのθに対応します。

このとき、相対位置と相対速度について、

${\bf x}_{ri} = (\frac{1}{2}l_{i}\sin{\theta_{i}},\frac{1}{2}l_{i}\cos \theta_{i})\\ {\bf v}_{ri} = \frac{\partial}{\partial t}{\bf x}_{ri} = (\frac{1}{2}l_{i}\dot{\theta}\cos{\theta_{i}},-\frac{1}{2}l_{i}\dot{\theta}\sin{\theta_{i}})$

また、絶対位置 ${\bf x}_{i}$ について、

${\bf x}_{0}={\bf x}_{r0}\\ {\bf x}_{1}=2{\bf x}_{r0}+{\bf x}_{r1}\\ {\bf x}_{2}=2{\bf x}_{r0}+2{\bf x}_{r1}+{\bf x}_{r2}$

速度 ${\bf v}_{i}$ も上式を微分するだけです。欲しいのは $\frac{1}{2}m_{i}|{\bf v}_{i}|^{2}=\frac{1}{2}m_{i}{\bf v}\cdot{\bf v}$ なので、展開して内積を計算すればそれぞれ導出できます。

回転のエネルギーについては $\dot{\theta}$ と慣性モーメントから計算できます。長さlの棒の慣性モーメントIは $I=\frac{1}{12}ml^2$ なので、回転のエネルギーについては $\frac{1}{24}ml^2\dot{\theta}^2$ と求まります。

なんやかんやで計算すると以下のようになります。

$T=(\frac{1}{6}m_{0}+\frac{1}{2}m_{1}+\frac{1}{2}m_{2})l_0^2\dot{\theta_0}^2\\ +(\frac{1}{6}m_{1}+\frac{1}{2}m_{2})l_1^2\dot{\theta_1}^2\\ +(\frac{1}{6}m_{2})l_2^2\dot{\theta_2}^2\\ +(\frac{1}{2}m_{1}+m_{2})l_0l_1\dot{\theta_0}\dot{\theta_1}\cos(\theta_0-\theta_1)\\ +(\frac{1}{2}m_{2})l_0l_2\dot{\theta_0}\dot{\theta_2}\cos(\theta_0-\theta_2)\\ +(\frac{1}{2}m_{2})l_0l_2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2)$

汚いようですがθに関しては対称的な形をしています。

同様に位置エネルギーに関しても

$U=(\frac{1}{2}m_0+m_1+m_2)gl_0\cos\theta_0+(\frac{1}{2}m_1+m_2)gl_1\cos\theta_1+(\frac{1}{2}m_2)gl_2\cos\theta_2$

ラグランジアンLは先に書いたようにL=T-Uとして求まります。あとはラグランジュの運動方程式

$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta_i}}-\frac{\partial L}{\partial \theta_i}=0$

を解けば $\ddot{\theta_i}$ について求まるという算段です。ただし、上の式を解こうとすると $\ddot{\theta_i}$ についての3元連立方程式になるので、微分方程式を解く前に適当な方法で変数を独立させておく必要があります。3変数なら一応解析的にも書き下すことができますが、大変面倒なので今回はガウス・ジョルダンの消去法を使いました。

ラグランジュ形式は便利ですが、最後に出てくる式がいまいち意味付けしにくいという不思議なところがあります。もっとも世の中の物理というのは別に人間にわかってもらうために運動しているわけではないのですが。
そもそもラグランジアンL自体なんの値なのかと考えだすと悩みどころはどんどん増えていきます。教科書には「意味のある値ではない」ときっぱり書いてあったりします。清々しいですね。

まあそれでもハミルトン形式まで持っていけばハミルトニアンHに意味を持たせることができるのですが。

運動方程式が解けたら次は微分方程式の数値計算の出番です。

微分方程式を数値計算で解く方法は数多くあります。今回のプログラムでは3種類から選べるようにしています。すべて陽解法で、それぞれ1次、2次、4次の精度を持っています。

おそらく最もシンプルでわかりやすい方法がオイラー法で、離散化した点の間を直線で補間する方法です。直線なので1次の精度がある、と表現されます。
オイラー法は精度が悪いため一般には使われません。左上のエネルギーの値のぶれが大きいことからもわかると思います。

より高い次数まで精度を持たせた方法があり、例えば4次のルンゲ・クッタ法が有名です。今回は2次のホイン法も用意しました。エネルギーの値を見れば違いがわかると思います。

ただし、たとえ4次精度を持っていたとしても長期間のシミュレーションを行うと次第にエネルギーが変化していってしまいます。長期間のシミュレーションでこういう性質があると結構問題で、例えば太陽系のシミュレーションをしていたら計算誤差で惑星が吹っ飛んでいったとか、分子の運動をシミュレーションしていたら計算誤差で勝手に沸騰した、なんていうよくわからないことが起きてしまいます。こういう場合にはエネルギーが長期間で安定するよう工夫した数値計算を行う必要があります。シンプレクティック性を持つ計算などと呼ばれますが、いちいち書いていたらとんでもない量になるのでここでは割愛します。

私の研究分野でよく使われる速度ベルレ法などの計算手法はうまい具合にこのシンプレクティック性を持っており、2次精度の手法ですが、計算が簡便でしかも陽解法でありながら長期間でエネルギーが変化しにくいという便利な性質を持っています。

系としてはとっても単純でありながら、運動はとても不規則です。時には大きくゆったり回ったり、時には先の振り子だけぐるぐる回ったり、それらを予告なく切り替え続けます。

このような運動はカオスと呼ばれますが、単純な法則から複雑な運動が生まれる例は他にもたくさんあって、例えば前回の反応拡散系も支配している方程式は単純ながら模様が複雑になったりもします。

それなら逆に、世の中の一見複雑な現象も実は単純な法則によって成り立っているのではないか?と考えるものの見方があります。では複雑な運動から単純な法則は見つけることができるのか?振り子ではなく、例えば世の中の気候の変化、株価、人の行動のような複雑な現象は実は単純な法則から成り立っているのではないか?そうだとしたらその法則は何か?という疑問に答えるのが、「複雑系」と呼ばれる学問のひとつの目標なのではないか、と私は勝手に考えています。

まあ、私は複雑系の研究者ではないので、好き勝手に妄想しているところが大きいのですが。

ぞうさんの何でもノート

思いついたことを書いていきます

月別アーカイブ: 2013年4月

[C++11対応記念] clang/LLVMのインストールとGMPを使って円周率ベンチ。

[作ってみた] [プログラミング] 剛体3重振り子