はじめに

Eigenライブラリを使って数値計算している場合、Eigenの提供する各種ソルバー

を入力データを基に動的に切り替えたいことがあります。
Eigenの密行列・疎行列用ソルバーは共通のインターフェースを持つため、テンプレートを使って容易に動的ポリモーフィズムを実現できます。

実装

追記：GitHub: shohirose/qiita/eigen_solversにアップしました。

solver_interface.hpp

#include <Eigen/Core>

// 線型方程式Ax = bを解くソルバークラスのインターフェース
template <typename MatrixType>
struct solver_interface {
  virtual ~solver_interface() = default;
  virtual void decompose(const Eigen::Ref<const MatrixType>& A) = 0;
  virtual Eigen::VectorXd solve(const Eigen::Ref<const Eigen::VectorXd>& b) = 0;
};

dense_solver_interface.hpp

#include "solver_interface.hpp"

using dense_solver_interface = solver_interface<Eigen::MatrixXd>;

sparse_solver_interface.hpp

#include <Eigen/SparseCore>
#include "solver_interface.hpp"

struct sparse_solver_interface : solver_interface<Eigen::SparseMatrix<double>> {
  using SparseMatrix = Eigen::SparseMatrix<double>;
  virtual ~sparse_solver_interface() override = default;
  virtual void analyze_pattern(const Eigen::Ref<const SparseMatrix>& A) = 0;
  virtual void factorize(const Eigen::Ref<const SparseMatrix>& A) = 0;
  virtual bool fail() const = 0;
};

dense_solver.hpp

#include <memory>
#include <string>
#include "dense_solver_interface.hpp"

template <typename SolverType>
class dense_solver : public dense_solver_interface {
 public:
  ~dense_solver() override = default;
  void decompose(const Eigen::Ref<const Eigen::MatrixXd>& A) override {
    solver_.compute(A);
  }
  Eigen::VectorXd solve(const Eigen::Ref<const Eigen::VectorXd>& b) override {
    return solver_.solve(b);
  }
 private:
  SolverType solver_;
};

std::unique_ptr<dense_solver_interface> make_dense_solver(const std::string& type);

sparse_solver.hpp

#include <memory>
#include <string>
#include "sparse_solver_interface.hpp"

template <SolverType>
class sparse_solver : public sparse_solver_interface {
 public:
  using SparseMatrix = Eigen::SparseMatrix<double>;
  ~sparse_solver() override = default;
  void decompose(const Eigen::Ref<const SparseMatrix>& A) override {
    solver_.compute(A);
  }
  void analyze_pattern(const Eigen::Ref<const SparseMatrix>& A) override {
    solver_.analyzePattern(A);
  }
  void factorize(const Eigen::Ref<const SparseMatrix>& A) override {
    solver_.factorize(A);
  }
  Eigen::VectorXd solve(const Eigen::Ref<const Eigen::VectorXd>& b) override {
    return solver_.solve(b);
  }
  bool fail() const override { solver_.info() != Eigen::Success; }
 private:
  SolverType solver_;
};

std::unique_ptr<sparse_solver_interface>
    make_sparse_solver(const std::string& type, const std::string& preconditioner = "");

dense_solver.cpp

#include "dense_solver.hpp"
#include <Eigen/Dense>

std::unique_ptr<dense_solver_interface> make_dense_solver(const std::string& type) {
  using namespace Eigen;
  using std::make_unique;

  if (type == "PartialPivLU")
    return make_unique<dense_solver<PartialPivLU<MatrixXd>>>();
  else if (type == "FullPivLU")
    return make_unique<dense_solver<FullPivLU<MatrixXd>>>();
  // ... 以下同様
  else {
    // エラー処理
  }
}

sparse_solver.cpp

#include "sparse_solver.hpp"
#include <Eigen/Sparse>

std::unique_ptr<sparse_solver_interface>
    make_sparse_solver(const std::string& type, const std::string& preconditioner) {
  using namespace Eigen;
  using std::make_unique;

  if (type == "BiCGSTAB") {
    if (preconditioner == "SimplicialCholesky") {
      return make_unique<sparse_solver<BiCGSTAB<
          SparseMatrix<double>, SimplicialCholesky<SparseMatrix<double>>>>>();
    } else if ( // 以下同様
       // ...
    } else {
       // エラー処理
    }
  } else if (type == "SparseLU")
    return make_unique<sparse_solver<SparseLU<SparseMatrix<double>>>>();
  else if (type == "SparseQR")
    // 同様
  else {
    // エラー処理
  }
}

使い方

密行列の場合

Eigen::MatrixXd A = // ...
Eigen::VectorXd b = // ...
auto solver = make_dense_solver("PartialPivLU");
solver->decompose(A);
const auto x = solver->solve(b);

疎行列の場合

Eigen::SparseMatrix<double> A = // ...
Eigen::VectorXd b = // ...
auto solver = make_sparse_solver("BiCGSTAB", "SimplicialCholesky");
solver->decompose(A);
// decompose()はanalyze_pattern()+factorize()と等価：
// solver->analyze_pattern(A);
// solver->factorize(A);
const auto x = solver->solve(b);

Eigenソルバーを動的に選択する

はじめに

実装

使い方