有限要素法(FEM)の実装におけるC++行列計算ライブラリEigenの使い方

最終更新 May 11, 2022 5 分で読める

C++行列計算ライブラリEigen

Eigenの使い方を調べている中で、日本語で大事な情報がないことが多かったので、特に疎行列に関して有限要素法(FEM)の実装で必要となる使い方と注意点をまとめた。

Eigenとは

C++で記述された線形代数ソルバーのオープンソースライブラリ
有名なGoogleのTensorFlowを含め、科学・ロボット・CG分野などで幅広く活用されている
テンプレートライブラリ、つまりヘッダファイルだけからなるので、 Eigen HomePageからダウンロードしてきたヘッダファイルをソースからインクルードするだけで利用可能.

// includeの例
#include <Eigen/Core>
#include <Eigen/Sparse>

参考となる資料

公式チュートリアルと公式クイックリファレンスは英語だが、一番情報が正確かつ網羅的なので困ったらこれらで調べる
日本語でもQiitaなどのネット記事がいくつかある

疎行列(Sparse Matrix)とは

行列の成分の大部分がゼロである行列のこと。 wikipediaの疎行列のページ

Eigenの疎行列の使い方

Eigen 3.3.7(2020/7/23の最新バージョン)をもとに、FEMの計算をEigenの疎行列で実装するときに必要な情報を中心に紹介する。

略称

Eigenの型が長いので、以下のようにtypedefする。

typedef Eigen::SparseMatrix<double, Eigen::RowMajor, int64_t> SpMat; // 疎行列の型
typedef Eigen::SparseVector<double> SpVec; // 疎ベクトルの型
typedef Eigen::MatrixXd Matrix; // 密行列の型
typedef Eigen::VectorXd Vector; // 密ベクトルの型

ここで注意すべき点は疎行列の格納の仕方を表すRowMajorである。デフォルトはColMajorであるが、ここではRowMajorとして記述していく (連立一次方程式の求解においてRowMajorのみがOpenMPに対応しているため)。

疎行列の宣言・サイズの指定・メモリの確保

int num_row = 100; // 行数
int num_col = 100; // 列数
int num_nonzero = 100; // 非ゼロ要素の数

SpMat spmat1(num_row, num_col); // 宣言+サイズの指定
spmat1.reserve(num_nonzero); // メモリの確保

SpMat spmat2; // 宣言
spmat2.resize(num_row, num_col); // サイズの指定
spmat2.reserve(num_nonzero); // メモリの確保

FEMにおいて非ゼロ要素の数はメッシュのコネクティビティの情報から事前に計算をしておく必要がある。

疎行列への値の代入

FEMにおいて要素行列から全体行列に足しこむ時に疎行列への値の代入が必要になる。ここでは2通りの方法を紹介する。

1. 直接代入する方法

spmat1.insert(0,0) = 1; // insert は計算速度は速いが既に要素が代入されている場合には使えない
spmat1.insert(0,0) = 2; // つまりここは更新されない
std::cout<<spmat1.coeffRef(0,0)<<std::endl; // この結果は1

spmat1.coeffRef(0,0) = 2; //  coeffRef は既に要素が代入されていても使えるが遅い
std::cout<<spmat1.coeffRef(0,0)<<std::endl; // この結果は2
spmat1.coeffRef(0,0) += 1; // coeffRef は値の加算ができるが計算速度は遅い
std::cout<<spmat1.coeffRef(0,0)<<std::endl; // この結果は3

**insertとcoeffRefの違いには注意が必要。 coeffRefは遅いので使うのはあまり推奨しない。** FEMにおいて全体行列を作成する場合には、事前に要素行列の和を計算してからinsertで代入する方法が速いと思われる。以下のようなデータ構造で要素行列の和を保持してから全体行列に足しこむのが1つの方法。

std::vector<std::vector<double>> idx; // idx[i][j] := i行目のj番目の非ゼロ要素の列番号
std::vector<std::vector<double>> adj_matrix; // adj_matrix[i][j] := i行目のj番目の非ゼロ要素の値

2. Tripletを使う方法

Tripletは(行番号, 列番号, 値)の組で管理するデータ構造であり、Eigenで定義されている、

typedef Eigen::Triplet<double> Triplet; // 略称
std::vector<Triplet> trivec; // Tripletを格納するvector
// Tripletを挿入していく.
trivec.push_back( Triplet(0,0,0.1) ); // 行番号:0, 列番号:0, 値:0.1
trivec.push_back( Triplet(1,1,0.2) ); // 行番号:1, 列番号:1, 値:0.2
// Tripletをまとめて疎行列に変換する
spmat.setFromTriplets(trivec.begin(), trivec.end());

Tripletを使う方法も計算速度は速いと思われるので、使いやすい方法で代入する方がよい。

行列ベクトル演算

行列ベクトル演算は、加減乗除が +, -, * , / で記述できるので簡潔に書ける。 MatrixXd, SparseMatrix, VectorXd, SparseVectorが混ざっていても問題なく計算できる。 FEMでは右辺項ベクトルの計算で疎行列ベクトル積が必要になる。

SpVec spvec(num_col); // 疎ベクトル
Vector vec1(num_col); // 密ベクトル
Matrix mat1(num_col, num_row); // 密行列
double coef = 1; // スカラー

// 値の代入は省略

SpMat spmat3 = spmat1 * spvec;   // 疎行列-疎ベクトル積
SpMat spmat4 = spmat1 + spmat2;  // 疎行列-疎行列和
SpMat spmat5 = spmat1 * coef;    // 疎行列-スカラー積
Vector vec2  = spmat1 * vec1;    // 疎行列-密ベクトル積は密ベクトルになる
Matrix mat2  = spmat1 * mat1;    // 疎行列-密行列積は密行列になる

疎と密の演算結果は密になるので注意が必要。

疎行列の要素のfor文による探索

FEMでは基本境界条件を設定する時にfor文による探索が必要となる。 (注意：Eigenの疎行列がRowMajorかColMajorかでfor文の回し方が異なる。以下はRowMajorの場合)

for(int i=0;i<A.outerSize();++i){ // outerSize()はcolumn-major matrixの列数,row-major matrixの行数
  for(SpMat::InnerIterator it(A, i); it; ++it){
    int id_col = it.col(); // 列番号
    int id_row = it.row(); // 行番号. ここではiと同じ
    // 基本境界条件の設定
    if(i == it.col()){
      it.valueRef() = 1; // 値の更新
    }else{
      it.valueRef() = 0; // 値の更新
    }
  }
}

連立一次方程式の求解

疎行列に対応している解法は疎行列対応ソルバーのリストに記載されている。また、マルチスレッド対応のEigenの計算は、Row-Major Sparse Matrix FormatのBiCGSTABである。並列化する場合は以下の1行が必要。

Eigen::initParallel();

BiCGSTAB法のコード(並列化対応)

// BiCGSTAB法(反復解法)の例
bool eigenBiCGSTAB(SpMat& A, Vector& x, Vector& b){
  Eigen::BiCGSTAB<SpMat> solver;
  solver.setTolerance(1e-9); // 許容誤差の設定
  solver.setMaxIterations(10000); // 最大反復回数の設定. デフォルトでは列数の2倍の数値
  solver.compute(A);
  if(solver.info()!=Eigen::Success) {
    // decomposition failed
    std::cout<<"Decomposition failed"<<std::endl;
    return false;
  }
  x = solver.solve(b);
  if(solver.info()!=Eigen::Success) {
    // solving failed
    std::cout<<"Failed solver BiCGSTAB"<<std::endl;
    return false;
  }
  std::cout << "#nbThreads:       " << omp_get_max_threads() <<std::endl;
  std::cout << "#iterations:     " << solver.iterations() << std::endl;
  std::cout << "estimated error: " << solver.error()      << std::endl;
  return true;
}

LU分解のコード

// LU分解(直接法)の例
bool eigenLU(SpMat& A, Vector& x, Vector& b){
  Eigen::SparseLU<SpMat> solver;
  solver.compute(A);
  if(solver.info()!=Eigen::Success) {
    // decomposition failed
    std::cout<<"Decomposition failed"<<std::endl;
    return false;
  }
  x = solver.solve(b);
  if(solver.info()!=Eigen::Success) {
    // solving failed
    std::cout<<"Failed solver LU"<<std::endl;
    return false;
  }
  return true;
}