「水素原子に対するSchrödinger方程式を有限要素法で数値的に解いてみる（C++のソースコード付き）」で紹介されているC++/Eigenのコードを、（ページ著者が）Juliaに移植してみたら12倍くらい速くなったと言うツイートが流れてきて、気になってC++の方のソースコードを見てみたのですが、
を活かしておらず、この2点を利用すればもっと時間を短縮できるものでした*2。実際、そう大きくない変更で、手元の環境では25分の1の時間で計算が終わるようになりました。
追記(2021/02/15 01:42):ブログ主が、有限要素法による水素原子のエネルギー固有値が、小さい方から見てどこまで正確かを検証する部分を加筆しました（図1, 図2）。これによって行列の固有値がすべて必要になるので、以下の方法は適切な計算方法にならないです。
なお、すべての固有値を得る場合にJulia比較でC++/Eigenが遅いのは、どうもC++/Eigenが自動的に並列化してマルチコアを活かしてこないところにあるので、その辺を何らかの方法*3で解消すると高速化できます。
EigenはC++用の線形代数ライブラリなのですが、名前がEigenなのに疎行列（Sparse Matrix）用の固有値問題のソルバーが無かったりします。上述のコードでは、GeneralizedSelfAdjointEigenSolverと言う計算する行列式に、疎行列、対称、正定値と言った条件を使わない上に、行列の固有値をすべて用いようとするソルバーが用いられています。汎用ソルバーとしても遅いというような指摘もあるようですが、入力データの特性を活かすのは高速化のイロハです。
この手の数値演算アルゴリズムを実装しようとすると数学の理解やデバッグの手間隙に大変な労力を費やすことになりますが、幸い、手軽に使えるテンプレート・ライブラリでSpectraと言うのがあります。具体的には、一般化固有値問題の計算対象の2つの行列をEigen::SparseMatrixで持ち、SymGEigsSolverと言うソルバーと解きます。なお、使ってみると癖があるなと思う部分があって、特にコンストラクタの第4引数を小さく取りすぎると固有値が収束しなかったとエラーを出して来るところには注意が必要です。
差分は以下になります。コンパイル時にダウンロードして展開したSpectraへインクルード・パスを通さないとコンパイル・エラーになるので気をつけてください。
--- a/src/hydrogen_fem/hydrogen_fem.cpp
+++ b/src/hydrogen_fem/hydrogen_fem.cpp
@@ -12,19 +12,24 @@
#include <memory> // for std::unique_ptr
#include <boost/assert.hpp> // for BOOST_ASSERT
#include <Eigen/Eigenvalues> // for Eigen::GeneralizedSelfAdjointEigenSolver
+#include <Eigen/Sparse>
+#include <Spectra/SymGEigsSolver.h>
+#include <Spectra/MatOp/SparseGenMatProd.h>
+#include <Spectra/MatOp/SparseCholesky.h>
+#include <iostream>
namespace hydrogen_fem {
// #region コンストラクタ
Hydrogen_FEM::Hydrogen_FEM()
- : hg_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL)),
+ : hg_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL)),
length_(ELE_TOTAL),
mat_A_ele_(boost::extents[ELE_TOTAL][2][2]),
mat_B_ele_(boost::extents[ELE_TOTAL][2][2]),
nod_num_seg_(boost::extents[ELE_TOTAL][2]),
node_r_ele_(boost::extents[ELE_TOTAL][2]),
node_r_glo_(NODE_TOTAL),
- ug_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL))
+ ug_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL))
{
}
@@ -42,20 +47,33 @@ namespace hydrogen_fem {
// 全体行列を生成
make_global_matrix();// 一般化固有値問題を解く
- Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> es(hg_, ug_);
+ Spectra::SparseSymMatProd<double> Aop(hg_);
+ Spectra::SparseCholesky<double> Bop(ug_);
+
+ // 第3引数の1は要求する固有値の数，第4引数は収束速度（とメモリー利用量と行列操作の複雑さ）を定める値
+ Spectra::SymGEigsSolver<double, Spectra::SMALLEST_ALGE, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, Spectra::GEIGS_CHOLESKY> geigs(&Aop, &Bop, 1, NODE_TOTAL<10 ? NODE_TOTAL : 10 + (NODE_TOTAL - 10)/100);
+
+ geigs.init();
+ int nconv = geigs.compute();
+ if(Spectra::SUCCESSFUL != geigs.info()){
+ std::cerr << "geigs.info(): " << geigs.info() << std::endl;
+ exit(-1);
+ }
// エネルギー固有値Eを取得
- auto const e = es.eigenvalues()[0];
+ Eigen::VectorXcd evalues = geigs.eigenvalues();
+ auto const e = evalues[0].real(); // 計算結果は複素数なことに注意．なお、複数の固有値を求めた場合、なぜか小さい順に固有値が並んでいるわけではない.
// 固有ベクトル（波動関数）を取得
- phi_ = es.eigenvectors().col(0);
+ phi_ = geigs.eigenvectors(1).col(0).real(); // なぜか添字が1からなのと、計算結果は複素数なのに注意
// 固有ベクトル（波動関数）を規格化
normalize();
return e;
}
void Hydrogen_FEM::save_result() const
@@ -190,8 +208,8 @@ namespace hydrogen_fem {
for (auto e = 0; e < ELE_TOTAL; e++) {
for (auto i = 0; i < 2; i++) {
for (auto j = 0; j < 2; j++) {
- hg_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
- ug_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
+ hg_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
+ ug_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
}
}
}
--- a/src/hydrogen_fem/hydrogen_fem.h
+++ b/src/hydrogen_fem/hydrogen_fem.h
@@ -15,6 +15,9 @@
#include <Eigen/Core> // for Eigen::MatrixXd
#include <boost/multi_array.hpp> // for boost::multi_array
+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+
namespace hydrogen_fem {
//! A class.
/*!
@@ -144,7 +147,7 @@ namespace hydrogen_fem {
/*!
左辺の全体行列
*/
- Eigen::MatrixXd hg_;
+ Eigen::SparseMatrix<double> hg_;
//! A private member variable.
/*!
@@ -192,7 +195,7 @@ namespace hydrogen_fem {
/*!
右辺の全体行列
*/
- Eigen::MatrixXd ug_;
+ Eigen::SparseMatrix<double> ug_;
// #region 禁止されたコンストラクタ・メンバ関数
