Lanczos法をRokkoのインターフェースに組み込む

[t-sakashita]

tutorial/titpack2にあるLanczos法のコードを、Rokkoの疎行列ソルバに組み込む。

• medium
• large

目的

• Anasazi、SLEPcがなくても使える疎行列のルーチンがほしい。
• TITPACK2の計算結果と比較したい。

[t-sakashita]

早めに行う。

[t-sakashita]

お願いします。 @wistaria

[t-sakashita]

• MPIプロセス数=1
• 簡単のため、Lanzosベクトルを行列で格納するので、非効率。

556a114d693cf7bf6b9bf90b6e757edd43a2c96a

[t-sakashita]

関数化した。 6defd1a004a2b4077a72a9b451d2a4c7535c3856 5c8f3e96a28207a5bdd69fcf7de892992549e786

[t-sakashita]

ScaLAPACKの三重対角行列に対する固有値計算ルーチンを用いる。

• alpha, betaは、全プロセスが全要素を持つ必要があるため、MPI通信によりgatherを行う。
• 固有ベクトルは、ブロックサイクリック行列として得られる。求める固有ベクトル数は数個であろう。固有ベクトル数 <= 列方向ブロックサイズと制限する。

[t-sakashita]

内積の並列版が必要。

template<typename T, int SIZE>
T inner_product(const Eigen::Vector<T, SIZE>& v, const Eigen::Vector<T, SIZE>& w, MPI_Comm comm) {
  T local_sum = v.dot(w);

  T sum;
  MPI_AllReduce(&local_sum, &sum, 1, rokko::mpi_type<T>,
                MPI_SUM, comm);

  return sum;
}

template<typename T, int SIZE>
T norm(const Eigen::Vector<T, SIZE>& v, MPI_Comm comm) {
  T local_sum = v.squaredNorm();

  T sum;
  MPI_AllReduce(&local_sum, &sum, 1, rokko::mpi_type<T>,
                MPI_SUM, comm);

  return sum;
}

da0923098ce50255837de92bf7ed97e2a2f69aa5

[t-sakashita]

現状では、全MPIプロセスで、同じ計算を実行している。 少ないメモリで巨大な行列の固有ベクトルを求めるためには、以下を行うべき。

• 固有値と固有ベクトルを求めるときで、別々にLanczos法を実行すべき。
• LanczosベクトルをEigen行列uに保存しない。 originalは本番計算用ではなく、コードの簡単なため、今のところ、上記は実装の改良はしない。