Intel Math Kernel Libraryでの大規模疎行列の直接解法＜動機＞
そこらの教科書による独学で力学分野における数値計算を学び、プログラム実装を行う場合、悲しいことに、連立方程式のソルバの検討に一番時間を食う。
たいていの場合、学習者のレベルが上がるにつれて、使用するソルバは
直接法（ガウスの消去法・LU分解）
↓
定常反復法（ガウス・ザイデル法、SOR法）
↓
非定常反復法（CG法ベース）
という流れになるのではないだろうか。日本語の教科書にはそのように書かれているし。計算対象が大規模になるにつれ、直接法より反復法、それもCGベースの非定常反復法が持て囃されていく。
しかしながら、CGベースの方法は対角成分は他より相対的に小さい場合には、まー収束しない。ひどい場合は発散する。この手の問題はだいたい何かの連成問題で頻繁に生じる（複数構造物の接触、流体構造連成、電気機械連成など）。ただでさえ複雑な定式化、実装を抱えて疲弊しているのに。
そんな理由から、数値計算として十分研究レベルになりうるレベルの計算を行う場合、普通の非定常反復法のさらにその先がほしくなる。その一番手っ取り早い方法がLU分解ベースの直接解法である。古くはスカイライン法やウェーブフロント法などがあるが、並列化効率のよいマルチフロンタル法が提案・実装されはじめ、欧米発のライブラリが結構世に出ている。Intel math kernel libraryに付属しているPardisoも大規模疎行列の直接解法ソルバーである。ここではその使用方法を備忘録としてまとめておきたい。ま、マニュアルが相当充実してるので、その必要もないかもしれないけど…。
静的ライブラリのコンパイルエラーの一例
（自分の中では）大規模な数値解析用プログラムを作成しており、静的ライブラリを複数作ってコンパイルする必要にかられている。
下記のリンクの問題を知らず、数時間ずっと悩んでいた。リンクの順番にも依存関係あるのね。
静的ライブラリのリンク時にundefined referenceエラーが出る(gcc)
知識がいい加減すぎるので、ちゃんと下記を読んで学ぶ次第。
- 作者: Robert Mecklenburg,矢吹道郎(監訳),菊池彰
- 出版社/メーカー: オライリージャパン
- 発売日: 2005/12/01
- メディア: 大型本
- 購入: 4人 クリック: 115回
- この商品を含むブログ (34件) を見る
Rでよく使うfunctionまとめ
よく忘れるので、完全に自分のためのメモを作成する。
クリップボードからデータのコピー
x=read.delim("clipboard")
散布図のプロット
plot(x,xlim=c(),ylim=c(),col="black",pch=2)
ピアソンの相関係数の計算
cor.test(x$dataname,x$dataname2,method="p")
回帰直線の計算
lm1<-lm(dataname2~datanamae1,data=x)
直交多項式（二次）
lm1<-lm(dataname2~poly(datanamae1,degree=2),data=x)
回帰線の描画をしたい場合
predict.c<-predict(lm1)
plot(X,predict.c,type="l")
回帰直線の作画
abline(lm1,lwd=1,col="blue")
図の上書き
par(new=T)
dataの全要素をcontensについて昇順に入れ替え
data2<-data[order(data$contents),]
gitのcommitエディタをvimに変更
情けないことにこれを調べずいつも戸惑っていた。
hikm.hatenablog.com
忘れないようリンク残します。
CentOS5.4 or 5.7へのgcc-4.8をインストール
やむを得ない事情で2016年現在でもcentOS5を使用しないといけない場合、
gccのバージョンが古すぎるためにコードのコンパイルができなかったり、
ソフトウェアのインストールができなかったりということが多々ある。
対応方法だが、下記のサイトに従うことで簡単にインストールできた。
www.task-notes.com
devtoolset-3はエラーが出たため諦めたが、devtoolset-2については問題ない。
なお、上記サイトの最後に/etc/profile.d内に作成するファイル名は、
devtools.shとしてshell scriptであることを明示しないと自分の環境では動作しなかった。
とにかくこれでc++11が好き放題使用できる。
C++からfortran呼び出し方法
爆速かつ可読性のよいコード作成のため、表題の内容に手を出した。自分の専門性を見失っている気がする。
基本的には下記を参考にすればよい。
kazuki-nagasawa.hatenablog.com
で、C++の場合は、プロトタイプ宣言を下記のようにCでexternしてやる必要がある。
extern "C"{
extern void add_array_( int*, int * );
extern void multiple_matrix_( double *, int *, int * );
};
また、コンパイルのリンクに-lstdc++を設定する。これだけで使用可能。
配列は次元に関わらず先頭のアドレスを渡す。また、定数であっても必ずアドレス渡しを設定すること。
もちろん動的配列も上記同様で、vector型ならば&v[0]、boost/multiarray型ならばv.origin()でよい。
なお、Cとfortranでは配列のアドレスのとり方が違うので、そこは注意がいる。
d.hatena.ne.jp
行優先か、列優先かは面倒なところだが、それ以外は相当簡単。
