C言語からLAPACKを呼んで逆行列を求める(LU分解する)
数値解析の授業で逆行列は求めちゃダメって言われたけど、気にしない気にしない。
実行列の逆行列を求める
Fortranの DGETRF でLU分解した後、 DGETRI で逆行列を求める。
/* * hoge.c * SIZE*SIZE型の実行列の逆行列を計算 * 元の行列は * (2 3 ) * (1 0.5) */ #define SIZE 2 #include <stdio.h> int main(void){ int m = SIZE ; // 行のサイズ int n = SIZE ; // 列のサイズ int lda = SIZE ; // mと同じ値 double A[SIZE*SIZE] ; // m x n の行列成分。この行列の逆行列を求める。 A[0] = 2.0 ;A[2] = 3.0; A[1] = 1.0 ;A[3] = 0.5; int info ; // 計算が成功すれば0を返す int ipiv[SIZE] ; // 要素数はm,nのうち小さい方とする int lwork = SIZE ; // nと同じ値 double work[SIZE] ; // 要素数はlworkと同じ値 // LAPACKのdgetrfサブルーチンを呼んで、行列AをLU分解 // 引数は全て参照渡し dgetrf_( &m, &n, A, &lda, ipiv, &info); // LU分解後の行列から逆行列を求める // 逆行列は元の配列Aに入る dgetri_( &n, A, &lda, ipiv, work, &lwork, &info); printf("%+8.5lf %+8.5lf\n", A[0], A[2]); printf("%+8.5lf %+8.5lf\n", A[1], A[3]); }
コンパイルと実行結果
$ gcc hoge.c -llapack -lblas -lm $ a.out -0.25000 +1.50000 +0.50000 -1.00000
引数の詳細については、下手な解説を見るより dgetrf と dgetri を読むのが一番わかりやすい。
CとFortranの配列形式の違いについては、 CとFortranで行列の添字が異なる点への注意喚起 に書いた。
複素行列の逆行列を求める
Fortarnの
ZGETRF
でLU分解した後、
ZGETRI
で逆行列を求める。
実行列の場合との違いは、行列Aと配列workの型がdouble
からdouble _Complex
に変わったところ。
/* * hoge.c * SIZE*SIZE型の複素行列の逆行列を計算 * 元の行列は * (0.5+i 1.0+0.5i) * (2.0 1.0 ) */ #define SIZE 2 // 2*2型の行列 #include <stdio.h> #include <complex.h> int main(void){ int m = SIZE ; // 行のサイズ int n = SIZE ; // 列のサイズ int lda = SIZE ; // mと同じ値 double _Complex A[SIZE*SIZE] ; // m x n の行列成分。この行列の逆行列を求める。 A[0]= 0.5+I ; A[2]= 1.0+0.5*I; A[1]= 2.0 ; A[3]= 1.0; int info ; // 計算が成功すれば0を返す int ipiv[SIZE] ; // 要素数はm,nのうち小さい方とする int lwork = SIZE ; // nと同じ値 double _Complex work[SIZE] ; // 要素数はlworkと同じ値 // LAPACKのzgetrfサブルーチンを呼んで、行列AをLU分解 // 引数は全て参照渡し zgetrf_( &m, &n, A, &lda, ipiv, &info); // LU分解後の行列から逆行列を求める // 逆行列は元の配列Aに入る zgetri_( &n, A, &lda, ipiv, work, &lwork, &info); printf("%+8.5lf%+8.5lfI %+8.5lf%+8.5lfI\n", creal(A[0]), cimag(A[0]), creal(A[2]), cimag(A[2])); printf("%+8.5lf%+8.5lfI %+8.5lf%+8.5lfI\n", creal(A[1]), cimag(A[1]), creal(A[3]), cimag(A[3])); }
コンパイルと実行結果
$ gcc hoge.c -llapack -lblas -lm $ a.out -0.66667-0.00000I +0.66667+0.33333I +1.33333+0.00000I -0.33333-0.66667I
引数の詳細は zgetrf と zgetri を読むべし。 またCとFortranの配列形式の違いについては、 CとFortranで行列の添字が異なる点への注意喚起 を読むべし。
Linux Mint にVim7.4をインストール
正式なaptリポジトリのVimはVer7.3で残念な感じだけど、ちょっと頑張ればVer7.4を入れることができる。
$ sudo add-apt-repository ppa:fcwu-tw/ppa
$ sudo apt-get update
$ sudo apt-get install vim
これでインストールできる。簡単だった。バージョンの確認は
$ vim --version |head VIM - Vi IMproved 7.4 (2013 Aug 10, compiled Aug 12 2013 09:20:12) Modified by pkg-vim-maintainers@lists.alioth.debian.org Compiled by buildd@ ...以下はコンパイルオプションの羅列...
+luaは嬉しいけど、-python3は悲しい。
Gvimもインストールする場合は、更に以下のコマンドを実行
$ sudo apt-get install vim-gtk
Ver7.4が使いたくて今まで自力でコンパイルしてきたけど、こんなに簡単に入るのか。あーあ
参考:
http://linuxg.net/how-to-install-vim-7-4-on-ubuntu-13-10-13-04-12-04-linux-mint-16-15-13-and-debian-sid/
https://gist.github.com/usufu/6398966
仮想環境のLinux MintにVmware Tools入れる時にkernel headersが見つからない問題
VMware Player にゲストOSとして Linux Mint を入れて、そこに VMware Tools をインストールするという話。 VMware Toolsを入れれば、仮想環境Windowを全画面にできるし、クリップボードが ホストOSとイケイケになるし、ディレクトリを共有できるし、かなり便利だ。
インストールの手順としては、ゲストOSの Linux Mint を立ち上げて、上タブの
仮想マシン
-VMware Toolsのインストール
をクリックすれば、DVDがマウントされる。
そのDVDの中に入っている VMware-Tools-8.8.6-12345.tar.gz
のような圧縮ファイルを解凍すれば vmware-tools-distrib
というディレクトリが生成される。
その中に vmware-install.pl
というperlスクリプトが入っている。
このvmware-install.pl
をsudoで実行して、エンターキーを押して進んでいくと、途中で
kernel headersのpathがどうたらこうたらで怒られてストップする(ことがある)。
解決方法
まず uname -a
して3.8.0-25-generic
のような文字列を探す。
$ uname -a Linux auewe-virtual-machine 3.8.0-25-generic #37-Ubuntu SMP Thu Jun 6 20:47:30 UTC 2013 i686 i686 i686 GNU/Linux
/usr/src/3.8.0-25-generic
以下にシンボリックリンクを貼る
/usr/src/linux-headers-3.8.0-26-generic/include/generated/uapi/linux/version.h
という実体へのシンボリックリンクを
/usr/src/linux-headers-3.8.0-26-generic/include/linux/version.h
に作成する。
$ sudo ln -s \ /usr/src/linux-headers-3.8.0-26-generic/include/generated/uapi/linux/version.h \ /usr/src/linux-headers-3.8.0-26-generic/include/linux/version.h
最後に vmware-install.pl
を実行
$ sudo perl vmware-install.pl
メモ:
export LANG=en
は意味なかった。
Python3の計算結果をGnuplotでグラフ化(subprocess Popenを使う)
ディスプレイにグラフを表示させる
pythonからgnuplotを使う上で、これが最も簡単な方法。
from subprocess import call call( ["gnuplot", "-e", "plot sin(x); pause -1"])
qを押せばグラフは消える。
グラフをファイルに保存
gnuplotの設定は長くなるので、 gnuplotCommand という文字列にまとめた。
from subprocess import call gnuplotCommand =''' set terminal png; set output "sin.png"; plot sin(x); ''' call( [ "gnuplot", "-e", gnuplotCommand])
Pythonで計算したデータをプロット
これが案外ハマる。 UNIX系なら echo で計算結果を標準出力に流して PIPE で gnuplot に食わせれば良いけれど、 Windowsでは echo が使えない(エラーが出る)。 代わりに Gow に含まれる printf 関数を使った。
from subprocess import Popen, PIPE data = ''' 1 5 2 10 3 15 ''' gnuplotCommand = ''' plot '-' ; pause -1 ''' printData = Popen( [ 'printf', data], stdout=PIPE) Popen( [ 'gnuplot', '-e', gnuplotCommand], stdin=printData.stdout)
シェルスクリプトからPython3に移行する人のために ~標準入出力・ファイル管理編~
汚いシェルスクリプトを保守する作業に嫌気がさしてきたので、Python3に乗り換えた。 せっかくなので必要になった知識を整理しておく。 同じようにシェルスクリフトから乗り換える人の役に立つと思う。
とりあえず IO 関連とファイル名の操作が重要な気がしたので、そのへんから書き始める。
環境変数の取得
# Python import os home = os.environ.get('HOME', "")
カレントディレクトリのパスを取得
# Bash ... もともと $PWDに入っているが current_path=`pwd`
# Python import os current_path = os.getcwd() # スクリプトのあるディレクトリの絶対パスは script_path = os.path.abspath(os.path.dirname(__file__))
ファイル名とディレクトリ名の結合
シェルスクリプトではファイルやディレクトリを文字列として扱い、 適当に切ったり貼ったりしていたけれど、 Pythonでは os.path.join を使って行儀よく扱う必要がある。
# Bash tmp_dir="$HOME/tmp" tmp_file="${tmp_dir}/tmp.txt"
# Python import os home = os.environ.get('HOME', "") tmp_dir = os.path.join( home , "tmp" ) tmp_file = os.path.join( tmp_dir, "tmp.txt") # 諸事情でパスの区切り文字を変更したくなったら tmp_dir_unix = tmp_dir.replace(os.path.sep, '/') # /path/to/tmp_dir tmp_dir_win = tmp_dir.replace(os.path.sep, '\\') # path\to\tmpdir
ディレクトリの作成
「同名のディレクトリの存在を確認して、無ければ作成」という定形作業
# Bash if [ ! -d ${tmp_dir} ] ;then mkdir ${tmp_dir} fi
# Python import os if not os.path.exists(tmp_dir): os.makedirs(tmp_dir) # 再帰的作成
ディレクトリの削除
いくつか方法があるけれど、一番強力なやつは
# Bash rm -rf ${tmp_dir}
# Python import shutil shutil.rmtree(tmp_dir)
ファイルへの書きこみ
# Bash echo "hoge\npiyo" > ${tmp_file}
#Python import codecs f_out = codecs.open(tmp_file, "w", "utf-8") print("hoge\npiyo", file=f_out) f_out.close()
ファイルの読みこみ
# Bash cat ${tmp_file}
# Python import codecs f_in = codecs.open(tmp_file, "r", "utf-8") for line in f_in.readlines(): print(line) f_in.close()
外部プログラムの実行とパイプ処理
# Bash cat ${tmp_file} | wc -l
# Python from subprocess import Popen, PIPE cat = Popen(["cat", tmp_file], stdout=PIPE) wc = Popen(["wc" , "-l"] , stdin=cat.stdout, stdout=PIPE) output = wc.stdout for line in output: print(line)
Python3の文字列について。%記法とformat()記法の対応関係
Python3では '%d'とか'%+5.3lf'のような記法は推奨されず、廃止される。 代わりにformat()を使う必要がある。 このような変更がなされた理由については、 Python 3.0 Hacks第6回 Pythonicな文字列フォーマットforamat()メソッド が分かりやすい。 高尚な理由はどうでもよくて、とりあえず昔の%記法との対応関係表が欲しいが見つからなかった。 仕方ないから、よく使うformatについてC言語のprintf関数との対応表を作った。 極めて愚直な表なので、とにかくこれさえ見れば新旧の違いが理解できる。
結果表示 | C言語 printf(%記法) | python print(format記法) |
---|---|---|
hoge(改行なし) | printf("hoge") | print("hoge", end="") |
hoge(改行あり) | printf("hoge\n") | print("hoge") |
hoge | printf("%s\n", "hoge") | print("{0}".format( "hoge")) |
hoge | printf("%s\n", "hoge") | print("{0:s}".format( "hoge")) |
hoge5piyo | printf("hoge%dpiyo\n", 5) | print("hoge{0}piyo".format( 5)) |
3 a 7 | printf("%d %s %d\n",3,"a",7) | print("{0} {1} {2}".format(3,"a",7)) |
------- | ------- | ------- |
5 | printf("%d\n", 5) | print("{0:d}".format( 5)) |
+5 | printf("%+d\n", 5) | print("{0:+d}".format( 5)) |
+005 | printf("%+04d\n", 5) | print("{0:+04d}".format( 5)) |
------- | ------- | ------- |
12.3456 | printf("%lf\n", 12.3456) | print("{0:f}".format( 12.3456)) |
_12.34 | printf("%6.2lf\n", 12.3456) | print("{0:6.2f}".format( 12.3456)) |
+12.34 | printf("%+6.2lf\n", 12.3456) | print("{0:+6.2f}".format( 12.3456)) |
012.34 | printf("%06.2lf\n", 12.3456) | print("{0:06.2f}".format( 12.3456)) |
------- | ------- | ------- |
_1.23e+01 | printf("%9.2e\n", 12.3456) | print("{0:9.2e}".format( 12.3456)) |
+1.23e+01 | printf("%+9.2e\n", 12.3456) | print("{0:+9.2e}".format( 12.3456)) |
01.23e+01 | printf("%09.2e\n", 12.3456) | print("{0:09.2e}".format( 12.3456)) |
6.2f
は、小数点2桁だけ必ず表示しつつ、6文字分の表示領域に収める場合(ただし999以上なら6文字オーバー)9.2e
は、小数点2桁だけ必ず表示しつつ、指数表記で9文字分の表示領域に収める場合
これぐらい知っとけば、一応なんとかなる。
メモ
その他、C言語にはない新しい表記法については、 format関数による文字列フォーマット(新しい形式) が分かりやすい。
公式の解説は 7.1.3. 書式指定文字列の文法
カーリーブレイス '{' を表示するには Python>str.format()はカールが使えない
Pythonで並列計算(Poolとmapを使う)
pythonで並列計算をしたいとき、Pool
とmap
を使うのが簡単だった。
C言語で並列計算して配列に格納する のと同じような感覚。
mapを使うのがカッコイイ。Haskellみたい。
簡単な例
- twiceFuncという名の関数 f(x) = 2x に
- argList = [1,2,3] の各要素を並列に代入して、
- resultList = [2,4,6] を得るプログラム
from multiprocessing import Pool def twiceFunc(arg): return 2 * arg if __name__ == '__main__': argList = [ 1, 2, 3] p = Pool() resultList = p.map(twiceFunc, argList) # resultList ... [ 2, 4, 6] print(resultList[0], resultList[1], resultList[2])
並列化された関数の終了を通知するテク
上記のコードのままでは、全計算が修了するまで結果が表示されないので不便。 そこで上記のtwiceFuncを以下のように変更すれば、各並列計算の終了がエラー出力に通知される。
from datetime import datetime def parallelCalcFunc(arg): # arg に計算をほどこして answer を作成 finStr = datetime.now().strftime("%Y/%m/%d %H:%M:%S") finStr += " arg:" + str(arg) finStr += " ans:" + str(answer) print(finStr, file=sys.stderr) return answer
参考
Python で並列計算 (multiprocessing モジュール) | 複数の引数を取る関数を map() メソッドで並列に走らせる