C++で行列演算を手軽に行うためによく使われれる（と思う）Eigenというライブラリがあります。純粋にヘッダファイルだけから構成されていて、C++コンパイラさえあればプラットフォームを選びませんし、それなりに高速性も考慮されていますし、縦横サイズも自由ですし、複素数の行列も扱えます。Eigenと名前が付いてるだけあって固有値、固有ベクトルの計算もほぼメソッド一発の楽ちんさんです。
ただ、これを使っていて変数宣言の型指定にauto(型推論)を使ったら数日ハマったことがあったのでメモを残しておきます。検証用コードは以下のとおりです。手元のcygwin g++ 4.9.2とMac OS Xのg++ (g++ 4.2.1 / LLVM 6.1.0)でおかしな事になるのを確認しています。
#include <iostream>
#include <complex>
#include "Eigen/Core"
#include "Eigen/LU"
#include <typeinfo>
using namespace std;
using namespace Eigen;
int main(void) {
MatrixXd m0(2, 2);
m0 << 1, 2, 3, 4;
const auto m1 = m0 * m0.inverse() * m0;
cout << "m1:\n" << m1 << endl;
const MatrixXd m2 = m0 * m0.inverse() * m0;
cout << "m2:\n" << m2 << endl;
}
処理内容としては何のこたぁない、2×2行列
m_0=\begin{pmatrix}1 & 2\\ 3 & 4\end{pmatrix}
に対して、$m_1, m_2$として$m_0\cdot m_0^{-1} \cdot m_0$を計算しようとしているだけで、coutによる出力としては両方とも$m_0$そのものが出力される事を期待します。が、cygwinとMac OS Xのg++では、$m_1$の値が実行する度に異なり、当然ながら期待した値になりません。一方で$m_2$の値は$m_0$と一致します。当たり前の期待値です。
$m_1$と$m_2$の違いはと言えば変数宣言時に型を型推論に任せている($m_1$)か、明示的に行列型をしている($m_2$)のみです。
・・・何が起きてるんでしょう? 公式ドキュメントやコンパイル結果、実行結果などから推察できるのは以下のようなことです。
Eigenでは行列計算時に遅延評価を用いています。そして行列(MatrixXd)同士の積演算の戻り値は行列(MatrixXd)型として定義されておらず、中間形式であるGeneralProduct型として定義されています。
本来ならEigenでの演算結果は$m_2$のように型を明示して宣言された変数への値格納が期待されています。$m_2$のように計算の結果を明示的に行列型変数に格納する場合には、暗黙の型変換により中間形式型から行列型へ値が変換され、その過程で値が評価・確定されます。
中間形式は演算の結果最終的に行列型の変数に格納されることが期待されており、$m_1$のように行列型変数格納前に該当行が終了してしまうと値を保持しているメモリが開放されてしまうようです。これが$m_1$の値が実行する度に変化してしまう理由でしょう。
どうしてもautoを使いたい場合でも、
const auto m1 = (m0 * m0.inverse() * m0).eval();
のようにeval()を用いるなど、複数の回避手段はあります。ただ、Eigenの思想に従うとするとEigenの行列を格納する変数をautoで宣言するのは避けるべきなのかな〜と思います。
わかってしまえばそういう思想のライブラリもあるよなとは思うのですが、「C++11以降はとにかくauto使おうぜ」ってあまり考えずに使っちゃうとはまりますのでご注意を。