Spectraとは
Spectraとは
ARPACKをC++用に再デザインされたもの。詳細は公式HPを参照。以下のような長所がある。
依存はEigenのみ。
ヘッダー・オンリーのライブラリ。
直感的なAPI。クラスを行列演算オブジェクトとしてソルバーに渡せる！
ARPACKをC++用に再デザインされたもの。詳細は公式HPを参照。以下のような長所がある。
依存はEigenのみ。
ヘッダー・オンリーのライブラリ。
直感的なAPI。クラスを行列演算オブジェクトとしてソルバーに渡せる！
ここでは、解きたい固有値問題をA v = ¥lambda vとする。行列AのサイズをNとする。Spectraには、固有値問題ソルバーとして、一般行列A用にGenEigsSolver、実対称行列A用にSymEigsSolver、エルミート行列A用にHermEigsSolverなどのクラスが用意されている。大規模行列Aを想定しているため、行列Aを直接与えるのではなく、演算Avを指定する任意のクラスopを引数として渡す。すなわち、Eigen SolverクラスはSpectra::GenEigsSolver<OpType> eigensolver(op, nev, ncv);のようにコンストラクトされる。クラスopは、以下のようなpublicメンバ関数を持たなければいけない。
using Scalar = type; //行列Aの要素の型（type = float, double, long double, or std::complex of them）を指定する。
int rows() const; //行列Aの行数を返すconst関数。
int cols() const; //行列Aの列数を返すconst関数。
void perform_op(const double *x_in, double *y_out) const; //行列演算y_out = A * x_inを定義する関数。
また、nevは計算する固有値の個数で、1 <= nev <= N - 2を満たす必要がある。ncvは収束速度を制御する整数である。大きなncvでは、収束が早くなるが、消費メモリが大きくなり各ステップでの行列演算が多くなる。nev + 2 <= ncv <= Nである必要があり、ncv >= 2 * nev + 1が推奨されている。
行列Aを保持しているなら、前述の演算クラスを簡単に定義できるヘルパークラス（see DenseGenMatProd, DenseSymShiftSolve, etc.）も用意されている。たとえば、DenseGenMatProd<Scalar>を用いると、
Eigen::MatrixXcd A;
A = ...;
DenseGenMatProd<std::complex<double>> op(A);
とするだけで行列演算オブジェクトop(A)を定義できる。op(A)は前述のpublicメンバ関数を有していることが簡単に確認できるはずだ。なおScalarには行列Aの要素の型（float, double, long double, そしてそれぞれのstd::complex）を与える。もちろんDenseSymMatProd<Scalar>にはstd::complexは指定できない。DenseHermMatProd<Scalar>にはstd::complexのみ指定可能。Eigen::SparseMatrix<type>に対する演算オブジェクトを定義する場合は、SparseGenMatProd, SparseHermMatProd, SparseSymMatProd, etc.を用いる。
例1) 行列Aを保持している場合。実対称行列Eigen::MatrixXd Aの最大固有値を3個求める。
#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
#include <iostream>
using namespace std;
using namespace Spectra;
int main()
{
// ランダムな実対称行列Aを生成する。
Eigen::MatrixXd M = Eigen::MatrixXd::Random(10, 10);
Eigen::MatrixXd A = M + M.transpose(); //対称化
// 行列演算オブジェクトのコンストラクト。
DenseSymMatProd<double> op(A);
cout << "Size of A: " << op.rows() << " x " << op.cols() << endl; //opはpublicメンバ変数にint rows() constやint cols() constなどを持つ。
// 固有値問題ソルバーのコンストラクト。ここではnev = 3個の固有値を求める。
SymEigsSolver<DenseSymMatProd<double>> solver(op, 3, 6);
// ソルバーの初期化。および計算実行。
solver.init();
int nconv = solver.compute(SortRule::LargestAlge);
// 結果の取得。
if (solver.info() == CompInfo::Successful) {
Eigen::VectorXd evalues = solver.eigenvalues();
cout << "Eigenvalues found:\n" << evalues << endl;
}
return 0;
}
例2) 行列演算y = A * xをクラスとして与える場合。
#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
#include <iostream>
using namespace std;
using namespace Spectra;
// A = diag(1, 2, ..., 10)とする。ただし行列Aは顕には定義しない。
class MyDiagonalTen {
public:
using Scalar = double; //doubleへのエイリアス。
int rows() const {return 10;} //行列Aの行数を返す関数。
int cols() const {return 10;} //行列Aの列数を返す関数。
// y_out = M * x_in
void perform_op(const double *x_in, double *y_out) const { //演算x_out -> y_out
for(int i = 0; i < rows(); i++) {
y_out[i] = x_in[i] * (i + 1);
}
}
};
int main() {
MyDiagonalTen op;
SymEigsSolver<MyDiagonalTen> eigs(op, 3, 6);
eigs.init();
eigs.compute(SortRule::LargestAlge);
if (eigs.info() == CompInfo::Successful) {
Eigen::VectorXd evalues = eigs.eigenvalues();
cout << "Eigenvalues found:\n" << evalues << endl; //ご存じのとおり、A = diag(1, 2, ..., 10)の固有値は1, 2, ..., 10です。
}
return 0;
}