B級科学者もどきの憂鬱

とある理系になりきれない奴のつれづれなる活動記

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

行列計算ライブラリEigenで固有値の計算

お久しぶりです。もうすぐ年明けですね。
もう何ヶ月も全然書いてなかったので、
次の更新はTIPSの話にしたかったのですが、
記事にできるほど進んでません!
というわけで今回は全然違うテーマです。

最近は数値シミュレーションを色々とやっています。
その過程で、行列計算にEigenというC++用ライブラリを導入しました。
色々見てる限り、結構速そうで導入も楽で使いやすそうです。

行列計算だけでなく、普通の配列を扱う計算で使っても、
アセンブラレベルで並列化したりしてくれるので、
普通にC++で書くより速いそうです。
これに関してはまだ確かめてませんけど。

導入は、公式サイトから圧縮ファイルをダウンロードし、
解凍して出来たファイルをインクルードパスに追加するだけです。
ヘッダファイルだけで出来ているので、コンパイル等は不要です。

例として、実数行列の固有値を求めるコードを書いてみます。

#include <stdio.h>
#include <Eigen/Dense>  // 固有値計算ルーチンの入ったインクルードファイル

int main() {
  Eigen::Matrix<double, 3, 3> A;  // 実数行列
  A << 1,2,3,4,5,6,7,8,9;  // Aの行列要素を代入

  Eigen::EigenSolver< Eigen::Matrix<double, 3, 3> > s(A);
  std::cout << "固有値\n" << s.eigenvalues() << std::endl;
  std::cout << "固有ベクトル\n" << s.eigenvectors() << std::endl;

  return 0;
}

EigenSolverが固有値を求めるテンプレートクラスで、
テンプレート引数に、求めたい行列の型名を入れて使います。
行列型の型名が割と長くなるので、typedef等で短い名前に置き換えた方が、
間違いが少なくなると思います。

複素行列の固有値を求めたい場合は、Aの型名の、
Eigen::Matrix<double, 3, 3>
doublestd::complex<double>に置き換えて、
(もちろん、EigenSolverのテンプレート引数も一緒に置き換えて)
Eigen::Matrix<std::complex<double>, 3, 3>
とすれば良い、だけではありません

これだけだと、EigenSolverのコンパイルに失敗します。
少なくとも、私が使った環境ではエラーになります。
(Eigen3.2.0、Visual Studio Express 2012 for Windows Desktop)
私は最初、これにハマって数十分悩みました。
というかこれが、今回の記事を書いた動機です。

解決法は、EigenSolverの代わりに、ComplexEigenSolverを使うことです。
これでエラーは出なくなります。

さらに詳しく知りたい方は、公式リファレンスを参照して下さい。

コメント

コメントの投稿


管理者にだけ表示を許可する

トラックバック

トラックバック URL
http://scientistb.blog42.fc2.com/tb.php/156-852efb75
この記事にトラックバックする(FC2ブログユーザー)

FC2Ad

まとめ

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。