Shifted BiCG法

とあるシミュレーションの問題では,シフト方程式と呼ばれる線形方程式
(A+\sigma I)x=b,\quad \sigma \in \mathbb{C},\quad A \in \mathbb{C}^{n \times n},\quad x,b \in \mathbb{C}^{n}
という,係数行列Aの対角成分がシフト量σだけ変化した方程式を解く必要が出てくるそうだ.Aは大規模で疎なものを対象としており,σは複数個現れるそうだ.これを解く方法としては,σだけ変化した係数行列をつくり,それをKrylov部分空間法なんかで解いてあげれば良いのだが,σの数だけ係数行列を新たに生成するのはあまり効率がよろしくない.

そこで,シフト方程式の係数行列を直接生成せずに,Ax=b(シード方程式と呼ばれている)を解く際に,ついでにシフト方程式も解いてあげて,全体の計算量を節約しよう,という解法が生み出された.さまざまな解法が生み出されているが,その1つに,BiCG法を拡張することで得られるShifted BiCG法がある.

論文に擬似コードが載っていたので,matlabで実装することにした.コード例はこんな感じ.シフト量σは複数個扱えるようにしている(はず).

シフト方程式が先に収束したらそれの反復は止めるようにしようとしたので,なんかとても中途半端な感じになっている(配列をまとめて処理するような書き方をしてしまったので,該当する要素については計算しない,というのがめんどくさくて投げた).そのうち直す(直さない).

Strassenのアルゴリズム

行列の積を求めるアルゴリズムに、シュトラッセンアルゴリズムと言うのがあるそうだ。

シュトラッセンのアルゴリズム - Wikipedia
Strassen algorithm - Wikipedia
Part II: The Strassen algorithm in Python, Java and C++ · Martin Thoma

まあ、wikipediaによると単純な3重ループより計算量が少し小さくなるらしい。

実装自体はそんなに難しくなさそうなので、scalaで書いてみた。


main関数は特に意味はない。テストのときはこれをimportして使った。

部分行列をどこまで分割していくかで結構性能が変わる。上の記事では一番小さい部分行列をleafと呼んでいたのでそれにならってここでもそのように表記する。

ためしに1024*1024の零行列同士の積を、leafの値を変えて計算し、かかった時間をグラフにした。
ちなみに使用したマシンは、数年前は新しかったiMac

データはこちら。

leafのサイズ[2^n * 2^n] 時間[ミリ秒]
0 318759
1 75572
2 25543
3 11921
4 6901
5 4567
6 3469
7 2956
8 2722
9 2737
10 3222

このマシンでは、leafの値は256ぐらいがよさそうだ。

ということで、どれくらい速くなるのかを

  • 単純な3重ループ(手計算でやるのと同じ手順。ijkと呼んでおく。)
  • 内側の2重ループの順番を入れ替えた3重ループ(ikjと呼んでおく。)
  • シュトラッセン(leaf=256)

の3つで、計算時間を比較した。

で、得られた結果をグラフにまとめたのがこちら。赤がシュトラッセン、青がijk、緑がikjとなっている。

ちょっとわかりにくいので行列のサイズが512のときのデータを示す。

アルゴリズム 時間[ミリ秒]
ijk 932
ikj 377
Strassen 368

シュトラッセン、あまり速くなさそうに見える。
ループの順番変えるだけでこんなに性能良くなるのには驚いた。

シュトラッセンアルゴリズムは、行列を次に大きい2^nの正方行列に拡大して計算するので、平べったい感じのグラフになる。なので、行列のサイズが130とか520とかだと、とても効率がよろしくないことになる。

ikjの方が色々な行列に対応できるので、扱う行列が2^nに近いものでなければ、こっちの方がいい感じなのかもしれない。



……なーんか腑に落ちないので、試しに大きめの正方行列でikjとシュトラッセンの時間を計ってみた。

行列のサイズ[N * N] ikj[ミリ秒] シュトラッセン[ミリ秒]
1024 3354 3121
2048 25765 21154

何回か計ってみたが、大体こんな感じになった。
この表を見ると、シュトラッセンの方が速いことが分かる。
大きい行列かつサイズが2^nに近いものであれば、シュトラッセンアルゴリズムは威力を発揮しそうだ。

Matrix Market形式の行列を読み込む

Matrix Market形式の行列を読み込んで2次元配列に突っ込んでみよう、というもの。

Matrix Market形式のファイルは、こんな感じになっている。

%hoge
%-------------------------------------------------------------------------------
148 148 1527
1 1 -99.99999999999983
2 1 24.999999999999957
3 1 7.839625797872487e-29
……

%はコメント行なので無視。
コメント以外の先頭の行には行列のサイズ、非零要素の数が書かれており、
そのあとに行番号、列番号、要素の値…と続いていく。

なので、ファイルを行区切りで読み込んだリストから行頭に%のないものを取り出し、先頭の要素の情報を使って配列を生成し、各要素を突っ込んでいくと、出来上がりとなる。

WinGW版OCamlのocamloptをcygwinなしで実行……

ここでのOCamlは公式サイトからダウンロードしたインストーラを使ってぶち込んだ4.00.1
MinGWはMinGW32(たぶん)

ocamloptというのを使ってみよう。

'i686-w64-mingw32-as' は、内部コマンドまたは外部コマンド、
操作可能なプログラムまたはバッチ ファイルとして認識されていません。


MinGW版のocamloptはcygwinを使わないと動かないらしい。

cygwin入れてやってみるかぁ〜。

動くには動くが、MinGWcygwinを両方入れておくのは気持ちが悪い。

cygwinのみでまとめてみようとするも色々とめんどくせぇことになる。

じゃあMinGWだけでできるのかしら?なんだかcygwin使わなくても動かしている記事をちょっと見かけた気がする。

MinGW\bin\にas.exeがあったので、名前似てるから代わりにこいつ走らせてみればいいんじゃね?ということで、i686-w64-mingw32-as.batを作成してパスの通っているところにぶん投げた。中身はこれだけ。

as %*


もう一回ocamloptをやってみる。

'i686-w64-mingw32-gcc' は、内部コマンドまたは外部コマンド、
操作可能なプログラムまたはバッチ ファイルとして認識されていません。
** Fatal error: Cannot run i686-w64-mingw32-gcc -print-sysroot


先ほどと同様にi686-w64-mingw32-gcc.batを作成。これも中身はこれだけ。

gcc %*


もう一回やってみる。

** Fatal error: Cannot find file "libws2_32"


そうだライブラリのパスを通していなかった。ということでオプションでパスを指定する。1箇所指定しただけだとまた同様のエラーが出たので、最終的に以下のような感じになった。

ocamlopt -ccopt "-L MinGW\libのパス" -ccopt "-L MinGW\lib\gcc\mingw32\4.7.0のパス" hoge.ml


動いた!!!一応!!!

……でもこれでいいのかしら?
これでやって壊滅的なことになってら泣く。

PowerShellでsayコマンドっぽいのを

やってみたいなと、ふと思ったので調べてみたらこんなのを見つけた。

PowerShell を使用してコマンドラインでしゃべらせる方法(.NET版) - CX's PowerShell Diary - PowerShellグループ

これをそのまま使わせていただこう、ということで関数にちょこちょこっと修正してプロファイルに書き足した。

ライブラリ読み込むと、コンソールにいろいろ情報が出力されるので、適当な変数に代入して出力しないようにした。自動変数$argsの中身をそのまましゃべる。

.tex->.pdfにするのにコマンド2回たたくのが面倒だからまとめる

タイトルのごとく、コマンドを続けて実行する関数を作った。

ファイル名を受け取り、platexとdvipdfmxを実行してpdfファイルを生成する。
pdfファイルができたら、ファイルを既定のプログラムで開く。
もしすでに何かのプログラムで同じ名前のファイルを開いていたりすると、そのファイルに変更を加えたりすることができなくなる場合があるので、関数の最初のほうで、生成されるpdfファイルと同じ名前のファイルを開いているプロセスを殺してから各コマンドを実行する。
また、オプションで-nologを付けると、pdfファイルができた後に.log、.aux、.dviを削除する。

以下、問題点。
ファイルと結びついてるプロセスを見つけるのに、ウィンドウタイトルを使っているので、そこにファイル名が出ないようなアプリケーションを使うと、この関数では対応できない。もっといい方法はないのかしら……
まあ自分はpdfを見るのにAdobe Readerしか使ってないから問題ないと言えば問題ないのだが……

偶数・奇数を判定するコード

twitterのタイムラインでたまにこいつが流れてくる。
偶数/奇数の判定 - Visual Basic 締切済み| 【OKWAVE】

書いてみた。

これを実行するとcのソースコードができる。こんな感じ。

#include<stdio.h>
#include<stdlib.h>

int a(b){
    int res;
    switch(b){
        case 0:
            res = 1;
            break;
        case 1:
            res = 0;
            break;
        case 2:
            res = 1;
            break;
        case 3:
            res = 0;
 
中略

        case 99:
            res = 0;
            break;
        default:
            printf("error\n");
            exit(-1);
    }
    return res;
}

int main(void){
    int n;
    while(1){
        scanf("%d", &n);
        printf("%d\n", a(n));
    }
    return 0;
}

といった感じのテキトーなものになる。
偶数だったら1、奇数だったら0が返ってくる。
いくつからいくつまでの偶奇を判定するかはPythonのコードにある変数iを変えてあげるといい。

いやー、しょうもないっすね。