March 04, 2009
自作行列ライブラリのGSL化
僕は車輪の再発明は大の得意(笑)で、『自作行列ライブラリ』などというものをこさえて使ってきました。しかしこの度、計算が怪しい(固有値分解で、対象が対角行列に近いと収束しないケースがある)や、計算速度をあげたい、という欲求に駆られましたので、表面では自作のライブラリを生かしつつ、裏側でごにょごにょして問題解決する、という方策をとってみることにしました。というと格好良く聞こえるかもしれませんが、実のところインターフェイスを変えてしまうと、今まで作成したプログラムのコードを全て変更しなければならず面倒くさいというのが本音です。
実はこの手の問題には今までチャレンジしたことがあります。というのもこの自作行列ライブラリをDSPに移植する際、高速化のために用意されたアセンブラルーチンに内部的に差し替えるということをしたことがあり、その際はC++言語の機能であるテンプレートの特殊化を使ってコードを差し替えました。
今回も同じ手法をとることにしますが、差し替え先のライブラリとしてGNU Scientific Library(GSL)を利用することにしました。行列演算のみならず、FFTや数値積分といった様々な算術アルゴリズムが用意されていること、OSを問わず使えること(Visual Studioでもコンパイル可能、Brian Gladmanさんが公開しているプロジェクトファイルを利用して構築可能)、Rubyへのポーティングがされていること(Ruby/GSL)、などが決め手になりました。行列についてはBlasのC版(CBlas)を基本にしているようですが、これを用途によって差し替えることによってさらにチューニングが可能であることも魅力です。
結果、できあがったのは元のmatrix.hに対するmatrix_gsl.h、complex.hに対するcomplex_gsl.h(チューニング途中であるため、今後も随時更新予定)です。テンプレートの特殊化によって、中身のアルゴリズムを入れ替えています。また、少しいいことがあってGSL自体はCのライブラリであるため、それ自体を扱おうとするとmalloc/freeをしっかりやる必要があるのですが、自作行列ライブラリについては参照カウンタを利用したフライウェイト・パターンになっているので、今回作成したライブラリを利用する限りはあまりメモリ管理を気にすることなく使えています。
アルゴリズムを実装することは理解も進んで僕にとっては楽しいのですが、一回つまり出すとなかなか前に進めないのが問題です。というわけで世の中の既存のツールをうまく活用できるようにもなりたいものです。
※2009/3/31追加、メモリが一部あふれていたので修正版に差し替えました。
