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

# はじめに

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

- [Dense linear problems and decompositions: Linear algebra and decompositions](http://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html)
- [Sparse linear algebra: Solving Sparse Linear Systems](http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html)

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

# 実装

追記：[GitHub: shohirose/qiita/eigen_solvers](https://github.com/shohirose/qiita/tree/master/shohirose/eigen_solvers)にアップしました。

```cpp: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;
};
```

```cpp:dense_solver_interface.hpp
#include "solver_interface.hpp"

using dense_solver_interface = solver_interface<Eigen::MatrixXd>;
```

```cpp: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;
};
```

```cpp: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);
```

```cpp: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 = "");
```

```cpp: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 {
    // エラー処理
  }
}
```

```cpp: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 {
    // エラー処理
  }
}
```

# 使い方

```cpp:密行列の場合
Eigen::MatrixXd A = // ...
Eigen::VectorXd b = // ...
auto solver = make_dense_solver("PartialPivLU");
solver->decompose(A);
const auto x = solver->solve(b);
```

```cpp:疎行列の場合
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);
```