疎行列での、任意のMPIコミュニケータ対応

MPIコミュニケータに保持するべきクラス

parallel_sparse_evで持つべき。
distributed_crs_matrixではない。

異なるMPIコミュニケータを複数を使いたい場合は、parallel_sparse_evを複数宣言するべき。

MPIコミュニケータの渡し方 https://github.com/trilinos/Trilinos_tutorial/wiki/EpetraLesson01

distributed_vectorは、ライブラリにもコミュニケータにも依存していない。

共通化したrokko::mapping_1dを復活させ、そのメンバ関数

num_local_rows

共通のmapping_1dクラスの作成

mapping_1dのメンバ関数とすべき：純粋仮想関数とする。

  int get_dim() const {
    return dim_;
  }
  int num_local_rows() const {
    return num_local_rows_;
  }
  int start_row() const {
    return start_row_;
  }
  int end_row() const {
    return end_row_;
  }

継承したクラスにおいて定義し、AnasaziやSLEPcの関数を呼び出す。その理由：

繰り返し使われる場合に、AnasaziまたはSLEPc用のオブジェクトを毎度、生成するのは時間がかかる。

distributes_crs_matrixのメンバ関数：

  int get_nnz() const {
    return matrix_->NumGlobalNonzeros();
  }

explicit distributed_crs_matrix(rokko::mapping_1d const& map)

呼び出し側の書き方：

  rokko::parallel_sparse_ev solver(library);
  auto map = solver.default_mapping(dim, comm);
  rokko::distributed_crs_matrix mat(map);

SLEPcでは、CRSを作った後でなければ、サイズの情報が入手できない。

MatGetOwnershipRange
SLEPcでも、mfreeでは、自分でサイズを指定する。

  MatCreateShell(PETSC_COMM_WORLD, mat->get_num_local_rows(), mat->get_num_local_rows(), mat->get_dim(), mat->get_dim(), mat, A);

そのため、rokko::slepc::mapping_1dでは、SLEPcの関数を用いない。

とりあえず、SLEPcでも、mapでstart_row, end_row, get_num_local_rowsを実装する。メインプログラムで使うのは、matのメンバ関数のままで、実装・検証を進めていく。

CRSでの行の分割は、決まっているので、共通化できるはず。

プロセス番号順に等分割
余りが生じた場合は、プロセス0から順に担当する。
検討：とびとびの行を担当させる等の拡張はあり得るか？

上記の計算を行うメンバ関数を持つデフォルト版のクラスが欲しい。これは、基底クラスrokko::mapping_1dとするか？その場合、純粋仮想ではなく、実装することになる。

各ソルバのmapping_1d

SLEPc向け：

  explicit mapping_1d(int dim, mpi_comm const& mpi_comm_in)
    : dim_(dim), mpi_comm_(mpi_comm_in), num_local_rows_(calculate_num_local_rows()) {
    rokko::mapping_1d::set_solver_name("slpec");
  }

solver_nameで処理を分ける。変換の失敗のNULLや例外の検出で対応できるならば、そうしたい。

共通のdistributed_crs_matrixのコンストラクタ

num_entries_per_rowは、行列によって異なる。そのため、mapping_1dではなく、distributed_crs_matrixのコンストラクタ引数として与えるべき。

引数の型は共通のrokko::mapping_1dでは判別できない。以下の引数を多重定義する。

rokko::anasazi::mapping_1d
rokko::slepc::mapping_1d

distributed_crs_matrixのコンストラクタでは、parallel_sparse_evのメンバsolver_nameをif文でチェックして、AnasaziかSLEPcかを判断する。

  explicit distributed_crs_matrix(rokko::mapping_1d const& map, int num_entries_per_row) : map_(static_cast<const rokko::anasazi::mapping_1d*>(&map)) {
    matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, map_->get_epetra_map(), num_entries_per_row));
    set_sizes();
  }

上記の方法ではダメ。

rokko::mapping_1dが基底型。
rokko::mapping_1d map = solver.default_mapping(dim, comm)で受け取る。
オブジェクトmapを作成する際に、派生型（rokko::anasazi::mapping_1dとrokko::slepc::mapping_1d）の情報が失われてしまう。
crs_impl_の型は、rokko::anasazi::mapping_1dまたはrokko::slepc::mapping_1d。
決め打ちで使いたい場合は、これらの派生型を直接作成し、distributed_crs_matrixのコンストラクタに渡せば良い。よって、mapping_1dもファクトリで生成するべき

各ソルバ向けのdistributed_crs_matrix

Anasazi向け

class distributed_crs_matrix : public rokko::detail::distributed_crs_matrix_base {
public:
  explicit distributed_crs_matrix(rokko::anasazi::mapping_1d const& map, int num_entries_per_row) : map_(dynamic_cast<const rokko::anasazi::mapping_1d*>(&map)) {
    matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, map_->get_epetra_map(), num_entries_per_row));
    set_sizes();
  }

private:
  const rokko::anasazi::mapping_1d* map_;

SLEPc向け

class distributed_crs_matrix : public rokko::detail::distributed_crs_matrix_base {
public:
  distributed_crs_matrix() {}
  ~distributed_crs_matrix() {}

  explicit distributed_crs_matrix(rokko::slepc::mapping_1d const& map, int num_entries_per_row) : map_(dynamic_cast<const rokko::slepc::mapping_1d*>(&map)) {
    initialize(map_->get_dim(), map_->get_dim(), num_entries_per_row);
  }

private:
  const rokko::slepc::mapping_1d* map_;

mfreeとの対応

CRSでも、get_local_offset()は必要か？

CRSと異なり、mfreeでは、ライブラリ依存の情報はユーザは全く意識せずに、呼び出すことができる。

  rokko::heisenberg_mfree mat(L, lattice);
  rokko::parameters info = solver.diagonalize(mat, params);

そのほか

mpi_vectorとの類似性

作業記録

  explicit distributed_crs_matrix(rokko::mapping_1d const& map) : map_(static_cast<const rokko::anasazi::mapping_1d*>(&map)) {

constを付ける：

  const mapping_1d* map_;

新インターフェスの完成後、以下を削除する。

  virtual rokko::detail::distributed_crs_matrix_base* create_distributed_crs_matrix(int row_dim,
										    int col_dim) = 0;
  virtual rokko::detail::distributed_crs_matrix_base* create_distributed_crs_matrix(int row_dim,
										  int col_dim, int num_entries_per_row) = 0;

diagonalize関数

現状では、SLEPc向けに作成したCRSをAnasaziで使ってもあらわなエラーとはならないはず。そこで、キャストが失敗した場合を検出する。

    A = reinterpret_cast<slepc::distributed_crs_matrix*>(mat.get_matrix())->get_matrix();

boost::variantの使用は検討しない。

ライブラリの判別（SLEPcか？　Anasaziか？）には、boost::variantは使用できない。その理由：

ファクトリでは、実行時に型が追加される。
boost::variantは静的

ファクトリの作成

ROKKO_REGISTER_PARALLEL_SPARSE_CRS(rokko::anasazi::distributed_crs_matrix, "anasazi", 40)

これで、if文による場合分けが不要となる。

mapping_1dもソルバによって異なるので、ファクトリにより生成するようにするべきかもしれない。

ファクトリにより生成すると、ポインタアクセスになるので、コストがかかりそう。

戻り値がAnasaziのオブジェクトなので、ファクトリ化できない。

  Teuchos::RCP<Epetra_CrsMatrix> get_matrix() {
    return matrix_;
  }

ファクトリ化できなくても、これまで通り、Anasazi用のラッパーにて以下のように取り出せる。

    problem_ = Teuchos::rcp(new eigenproblem_t(reinterpret_cast<anasazi::distributed_crs_matrix*>(mat.get_matrix())->get_matrix(), multivector));

AnasaziおよびSLEPcのラッパーで、rokko::distributed_crs_matrixをそれぞれの型に変換する関数をオーバーロードした。

  parameters diagonalize(rokko::distributed_crs_matrix& mat, rokko::parameters const& params) {
    return diagonalize(static_cast<rokko::slepc::distributed_crs_matrix&>(mat.get_matrix()), params);
  }

CRSのtraitsは要らない。

以下は、ダメだった。巡回参照を引き起こすため、リンクがうまくいかない。

template<typename SOLVER>
class ps_ev_wrapper : public ps_ev_base {
  using solver_type = SOLVER;
  using crs_type = rokko::crs_t<solver_type>;
public:
  ps_ev_wrapper() : solver_impl_() {}
  virtual ~ps_ev_wrapper() {}
  void initialize(int& argc, char**& argv) { solver_impl_.initialize(argc, argv); }
  void finalize() { solver_impl_.finalize(); }
  parameters diagonalize(rokko::distributed_crs_matrix& mat, rokko::parameters const& params) {
    return solver_impl_.diagonalize(static_cast<crs_type&>(mat.get_matrix()), params);
  }

メインプログラムの書き方

Anasazi専用

  rokko::anasazi::solver solver;
  rokko::anasazi::mapping_1d map(dim, rokko::mpi_comm{MPI_COMM_WORLD});
  rokko::anasazi::distributed_crs_matrix mat(map, 10);

  rokko::parameters info = solver.diagonalize(mat, params);

SLEPc専用

  rokko::slepc::solver solver;
  rokko::slepc::mapping_1d map(dim, rokko::mpi_comm{MPI_COMM_WORLD});
  rokko::slepc::distributed_crs_matrix mat(map, 10);

  rokko::parameters info = solver.diagonalize(mat, params);

  explicit distributed_crs_matrix(rokko::mapping_1d const& map, int num_entries_per_row) : map_(dynamic_cast<const rokko::slepc::mapping_1d*>(&map)) {
    initialize(map_->get_dim(), map_->get_dim(), num_entries_per_row);
  }
  explicit distributed_crs_matrix(rokko::slepc::mapping_1d const& map, int num_entries_per_row) : map_(dynamic_cast<const rokko::slepc::mapping_1d*>(&map)) {
    initialize(map_->get_dim(), map_->get_dim(), num_entries_per_row);
  }

solverクラスで、名前を返すsolver_name()を利用する。

結果的に、以下のように、基底クラスを共通化した。

rokko/mapping_1d.hpp

class ps_mapping_1d_wrapper : public ps_mapping_1d_base {

rokko/anasazi/mapping_1d.hpp

class mapping_1d : public ps_mapping_1d_base { // ファクトリを使用するので、この継承は不要である。ポインタをダウンキャストするために必要。

std::shared_ptrのキャストには、std::dynamic_pointer_castを使う。

std::shared_ptrから生ポインタを取り出す。

  explicit distributed_crs_matrix(rokko::mapping_1d const& map, int num_entries_per_row) : map_(dynamic_cast<rokko::slepc::mapping_1d*>(map.get_impl().get())) {

rokko::mapping_1dでファクトリを使用した。

  rokko::parallel_sparse_ev solver(library);
  rokko::mapping_1d map(dim, rokko::mpi_comm{MPI_COMM_WORLD}, solver);

以下は、rokko::parallel_sparse_evが巡回依存になるため、コンパイルできない。

  mapping_1d(int dim, mpi_comm const& mpi_comm_in) : mapping_1d(dim, mpi_comm_in, rokko::parallel_sparse_ev::default_solver()) {}

https://stackoverflow.com/questions/29141629/c-generic-factory-that-can-call-any-constructor/29154449

https://stackoverflow.com/questions/23012621/pure-virtual-and-stdshared-ptr

In file included from /Users/sakashitatatsuya/development/rokko/rokko/anasazi/core.hpp:19:
/Users/sakashitatatsuya/development/rokko/rokko/anasazi/mapping_1d.hpp:30:12: error: constructor for
      'rokko::anasazi::mapping_1d' must explicitly initialize the member 'map_' which does not have a default constructor
  explicit mapping_1d(int dim, mpi_comm const& mpi_comm_in) {
           ^
/Users/sakashitatatsuya/development/rokko/rokko/anasazi/mapping_1d.hpp:63:14: note: member is declared here
  Epetra_Map map_;

ポインタのキャストが失敗していないかをチェックする：

    const rokko::slepc::distributed_crs_matrix* ptr = dynamic_cast<const rokko::slepc::distributed_crs_matrix*>(mat.get_ptr()->get_impl());
    if (ptr == nullptr)
      throw std::invalid_argument("SLEPc's diagonalize() : Not SLEPc's distributed_crs_matrix is given.");

以下の関数にconstを付けて良いか？独自のスマートポインタが戻り値である。ポインタはコピーされるか？

  Teuchos::RCP<Epetra_CrsMatrix> get_matrix() const {
    return matrix_;
  }

Epetra_MpiCommとEpetra_Mapにはデフォルトコンストラクタがない。そのため、ファクトリにより、rokko::anasazi::mapping_1dやrokko::slepc::mapping_1dを生成する際に困るはず。

ep_comm_(new Epetra_MpiComm(mpi_comm_in.get_comm())), map_(new Epetra_Map(dim_, 0, *ep_comm_))

以下は不要：

  const rokko::anasazi::mapping_1d& get_map_obj() const { return *map_; }

ポインタから参照を作る。

    rokko::anasazi::distributed_crs_matrix* ptr = dynamic_cast<rokko::anasazi::distributed_crs_matrix*>(mat.get_ptr()->get_impl());
    if (ptr == nullptr)  std::cout << "nullptr" << std::endl;
    //rokko::anasazi::distributed_crs_matrix& m = *ptr;
    return diagonalize(*ptr, params);

rokko::mapping_1dを生成する方法：

  rokko::mapping_1d map = solver.default_mapping(dim, rokko::mpi_comm{MPI_COMM_WORLD});

  rokko::mapping_1d map(dim, rokko::mpi_comm{MPI_COMM_WORLD}, std::string("anasazi"));

以下は、不採用：

  rokko::mapping_1d map(dim, rokko::mpi_comm{MPI_COMM_WORLD}, solver);

済

solverのメンバ関数default_mappingを定義する。
mapping_1dの基底クラスでのサイズ・添字計算を汎用化。SLEPcで用いる。
rokko::mpi_comm mpi_comm_は、anasaziとslepcの継承クラスで定義する。
以下は取りやめ。
- Anasaziではoverride
- 基底クラスps_mapping_1d_baseは純粋仮想ではなく、num_local_rows_やdim_を持つ。
- num_local_rows_やdim_の初期化には、継承コンストラクタを用いる。（基底クラスでprotected属性がついていても、派生クラスでメンバ初期化を行うことはできない。）→コンストラクタに与える変数の順序を覚えなければならない。
デフォルトソルバを扱い。（空文字""を書く必要はない。）

  mapping_1d(int dim, mpi_comm const& mpi_comm_in) : mapping_1d(dim, mpi_comm_in, detail::ps_mapping_1d_factory::instance()->default_product_name()) {}

以下の継承をやめた。

ps_mapping_1d_base : public rokko::mpi_comm

has-aの関係になるように、基底クラスps_mapping_1d_baseのメンバとして持たせた。

rokko::mapping_1dのメンバとして、sovler_name_を持つ。 rokko::anasazi::mapping_1dやrokko::slepc::mapping_1dのメンバとして持つと、間接アクセスで呼び出せるように、get_solver_name()をこれらのクラス全てで定義するのが面倒。

rokko::distributed_crs_matrixでget_solver_name()を定義し、その中で、rokko::anasazi::mapping_1dを定義する。

solver.diagonalizeでソルバ名が一致するかをチェックする。

solverクラスのget_solver_name()
rokko::distributed_crs_matrixのget_solver_name() このチェックをすれば、dynamic_castがnullptrになるかをチェックしなくても良い。

  parameters diagonalize(rokko::distributed_crs_matrix& mat, rokko::parameters const& params) {
    assert( get_solver_name() == mat.get_solver_name() );

文字列を比較する。エラーメッセージを詳しくできる。

  parameters diagonalize(const rokko::distributed_crs_matrix& mat, rokko::parameters const& params) {
    if (mat.get_solver_name() != "slepc") {
      throw std::invalid_argument("SLEPc's diagonalize() : " + mat.get_solver_name() + "'s distributed_crs_matrix is given.");
    }
    return diagonalize(*static_cast<const rokko::slepc::distributed_crs_matrix*>(mat.get_ptr()->get_impl()), params);
  }

    A = const_cast<Mat*>(mat.get_matrix());

laplacian_crs_mpiでのチェック

以下は不要：

  //void init(int dim, mpi_comm const& mpi_comm_in) { map_impl_->init(dim, mpi_comm_in); }

廃止した。

anasaziのdistributed_crs_matrixのset_sizes()
AnasaziとSLEPcのdistributed_crs_matrixのprivate: int num_local_rows_, start_row_, end_row_;を廃止するか？

AnasaziとSLEPcのmapping_1dから、num_local_rows_を削除した。その場で計算するか、ライブラリの関数を呼び出す。

エイリアスを止める。rokko::detail::distributed_crs_matrix_base

root_proc = 0をconstexprで定義

以下の仮引数のポインタをconstにしても、呼び出し元には制約は課されないはず。

  void insert(int row, int col_size, int* cols, double* const values) {

https://stackoverflow.com/questions/49652710/overload-class-member-function-marked-const

引数int row_dim, int col_dimをstd::arrayとした。 bc2c1b1ffd35bf63cbb16f2222011fec3b318e6c

委譲コンストラクタを書いてみた。

  template<typename SOLVER>
  distributed_crs_matrix(std::tuple<int,int> const& dims, SOLVER& solver_in) : distributed_crs_matrix(rokko::to_array(dims), solver_in) {}

引数{100,100}を渡すと、std::tuple<int,int>とstd::array<int,2>のどちらのオーバーロードを使うかが曖昧なため、エラーとなる。そのため、Python専用ラッパーを書くことで対応する。

Pythonのwrap_communicatorを、rokko::mpi_commの継承で書き直した。

関数create_distributed_crs_matrixはconstとした。staticにはできない。

mfreeに対して、Fortran用の1始まりのstart_row、end_rowのgeneric nameを定義した。

やるべきこと

メンバとして保持するポインタには、std::shared_ptrを用いる。
MPIコミュニケータを任意にする。
solverクラスにおいて、MPI_Commを保持

関数をrokko::anasazi::mapping_1d const*に帰着させる。

  void initialize(rokko::anasazi::mapping_1d const& map, int num_entries_per_row) {

distributed_mfreeで、mapping_1dを使う。任意のコミュニケータは対応済

Epetra_CrsMatrixは長方行列に対応している。

NumEntriesPerRowをcol_dimとするのは正しいか？

mapping_1dは１次元なので、col_dimは扱えない。

insert関数は、INSERT_VALUESに対応させる。
ADD_VALUESに対応する関数を作成する。

出力

ソルバの進捗情報はデフォルトoff

「ポインタの参照渡し」が使えるか？

CRSの行を取り出すインターフェース関数を取り出す。

py::initでラムダ関数が使える。 https://stackoverflow.com/questions/52657173/sharing-an-mpi-communicator-using-pybind11

    .def(py::init([](int dim, pybind11::handle const& comm_handle, std::string const& solver_name)  {
          const int rc = import_mpi4py();
          assert(rc==0);
          wrap_communicator wrap_comm(*PyMPIComm_Get(comm_handle.ptr()));
          MPI_Comm comm = wrap_comm.get_comm();
          return new mapping_1d(dim, rokko::mpi_comm{comm}, solver_name);
        }))

以下を追加する：

  virtual MPI_Comm get_comm() const { return MPI_COMM_WORLD; }

rokko::mpi_commの方が良いか？他のコミュニケータを使いたい場合は、overrideする。 overrideを付けると、pure virtualも含めて）オーバーライドする全てのメンバ関数に、overrideをつけないと警告が出る。

mapping_1dのC/Fortranバインディングを作成する。

Pythonバインディングにおいて、mat.start_rowを、map.start_rowに書き換える。他の言語でも同様。

出力のmapが最後の引数となる。
コンストラクタでは、出力であるmatを先頭に書く。

  call rokko_default_mapping(solver, dim, mpi_comm_world, map)
  call rokko_construct(mat, map, L+1)

任意のMPIコミュニケータでのテスト

Pythonインターフェースでは、std::functionからmfreeを作成した。 C++インターフェースにおいても、std::fucntionからmfreeを作成できるようにした。

  distributed_mfree(std::function<const double*, double*> const& multiply, int dim, MPI_Comm comm)

Eigen::Vector版も作成した。
std::vector版は作成できない。（std::vectorにポインタを割りつけることができないため。）
引数として渡すと、std::fucntion型のメンバとして保持し、継承したmultiply関数の中で呼び出すことにした。（この呼び出しの実行効率は良いか？　インライン化はされるか？）
クラスとして継承することができない。
ユーザがmfreeを定義した方が実行効率は良いだろう。
dim, commはskel::mapping_1dで束ねて作成するようにした。ただし、そのように作成したmapping_1dは、multiplyから参照できないため、勿体無い。そのようなmapping_1dはskel::mapping_1dに限らなくても良い。利点はなく、その必要もないが、以下も使えるようにできる：
rokko::anasazi::mapping_1d
rokko::slepc::mapping_1d
rokko::mapping_1d（パラメータにより選べる。）

  const int dim = 1 << L;
  rokko::skel::mapping_1d map{dim, rokko::mpi_comm{MPI_COMM_WORLD}};
  std::vector<double> buffer(map.get_num_local_rows());
  std::function<void(const double *const x, double *const)> multiply = [&, L, lattice, buffer](const double *const x, double *const y) {
    rokko::heisenberg_hamiltonian::multiply(MPI_COMM_WORLD, L, lattice, x, y, const_cast<double*>(buffer.data()));
  };

キャプチャしたbufferのポインタには、const_castが必要だった。

  rokko::distributed_mfree_adaptor mat(multiply, rokko::mapping_1d(dim, rokko::mpi_comm{MPI_COMM_WORLD}));

https://stackoverflow.com/questions/4665117/c-virtual-function-return-type

不要なxとyを除く：

std::function<void(const double *const x, double *const y)

キャストすべき型を、distributed_mfreeからdistributed_mfree_cへ変更する。

struct rokko_parameters rokko_parallel_sparse_ev_diagonalize_distributed_mfree(struct rokko_parallel_sparse_ev solver,
									       struct rokko_distributed_mfree mat,
									       struct rokko_parameters params) {
  struct rokko_parameters params_out;
  rokko_parameters_construct(&params_out);
  std::cout << "MFFFFREEEE" << std::endl;
  *static_cast<rokko::parameters*>(params_out.ptr) = static_cast<rokko::parallel_sparse_ev*>(solver.ptr)->diagonalize(*static_cast<distributed_mfree_c*>(mat.ptr),

distributed_mfree_cやdistributed_mfree_fをヘッダファイルに収める。

検討

std::shared_ptrからget()関数で生ポインタを取り出すのは、避けることができた。

継承することの是非を検討する：

mapping_1d : public detail::ps_mapping_1d_base
class distributed_crs_matrix : public rokko::detail::distributed_crs_matrix_base

RokkoInterfaceSparseSolversArbitratyMPICommnunicator - t-sakashita/rokko GitHub Wiki

疎行列での、任意のMPIコミュニケータ対応

MPIコミュニケータに保持するべきクラス

共通のmapping_1dクラスの作成

各ソルバのmapping_1d

共通のdistributed_crs_matrixのコンストラクタ

各ソルバ向けのdistributed_crs_matrix

mfreeとの対応

そのほか

作業記録

diagonalize関数

boost::variantの使用は検討しない。

ファクトリの作成

メインプログラムの書き方

済

やるべきこと

検討

⚠️ GitHub.com Fallback ⚠️

RokkoInterfaceSparseSolversArbitratyMPICommnunicator - t-sakashita/rokko GitHub Wiki

疎行列での、任意のMPIコミュニケータ対応

MPIコミュニケータに保持するべきクラス

共通のmapping_1dクラスの作成

各ソルバのmapping_1d

共通のdistributed_crs_matrixのコンストラクタ

各ソルバ向けのdistributed_crs_matrix

mfreeとの対応

そのほか

作業記録

diagonalize関数

boost::variantの使用は検討しない。

ファクトリの作成

メインプログラムの書き方

済

やるべきこと

検討

⚠️ **GitHub.com Fallback** ⚠️

⚠️ GitHub.com Fallback ⚠️