Spectra

Spectraとは

• Spectraとは、C++用の大規模固有値問題ソルバーライブラリ(https://spectralib.org/)
• 構造解析コードで線形座屈解析を実装する時に、一般化固有値問題を解く必要があり、C++で一般化固有値問題を高速に解く疎行列対応のソルバーを探した結果Spectraが候補として見つかった。（Eigenにも固有値問題を解くソルバーはあるが遅い）
• 使い方について調べても日本語の情報が少なく、使えるまでに時間がかかったため、メモとして残す。

使い方

• 自作構造解析ソルバー (https://github.com/kwkm0429/fem_structure) に実装する際に、 参考文献2を参考に実装したが、エラーが出てしまった。
• Spectra::SymGEigsSolverの引数に関するエラーメッセージだったため、引数の型を調べながら変更した結果エラーが出ずに正常に動作した。
• 以下に実装にハマったポイントを示す。
  • 元々EigenのColMajorで疎行列を扱っていたが、Spectraではエラーが出たため、RowMajorに変更した
  • Eigenの疎行列AとBをSpectra::SparseSymMatProdとSpectra::SparseCholeskyに変える
  • 一般化固有値問題を解くソルバーとして、Spectra::SymGEigsSolverを定義するが、この時に引数とする型の指定を公式のドキュメントを調べて試行錯誤した結果、以下の型指定で上手くいった。
• 以下の実装したSpectraSolverという関数は、対称行列A、正定値行列Bに対して、Aφ=λBφを解く関数である。（SpMatなどはEigenの疎行列を別の所でtypedefしているため、関数だけコピペではエラーが出る。）

bool SpectraSolver(SpMat& A, SpMat& B, Matrix& phi, std::vector<double>& lambda, Sim& sim){
  int i, j, id, size;
  // convert to Eigen::ColMajor because Spectra and Eigen has error for RowMajor
  Eigen::SparseMatrix<double> Ae, Be;
  Ae.resize(A.outerSize(), A.outerSize());
  Be.resize(B.outerSize(), B.outerSize());
  Ae.reserve(sim.num_nonzero);
  Be.reserve(sim.num_nonzero);
  for(i=0;i<A.outerSize();++i){
    for(SpMat::InnerIterator it(A, i); it; ++it){
      Ae.insert(it.col(), i) = it.valueRef();
    }
  }
  for(i=0;i<B.outerSize();++i){
    for(SpMat::InnerIterator it(B, i); it; ++it){
      Be.insert(it.col(), i) = it.valueRef();
    }
  }
  Spectra::SparseSymMatProd<double> Aop(Ae);
  Spectra::SparseCholesky<double> Bop(Be);
  Spectra::SymGEigsSolver<Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, Spectra::GEigsMode::Cholesky> geigs(Aop, Bop, sim.num_mode+sim.num_spc, sim.num_mode*2+sim.num_spc);
  geigs.init();
  int nconv = geigs.compute(Spectra::SortRule::LargestAlge);

  size = geigs.eigenvalues().size();
  id = 0;
  std::cout<<geigs.eigenvalues()<<std::endl;
  for(i=sim.num_spc;i<sim.num_spc+sim.num_mode;i++){
    lambda[id] = 1.0/geigs.eigenvalues()[i];
    phi.col(id) = geigs.eigenvectors().col(i);
    std::cout << id<<") Eigen value, lambda = " << lambda[id] << std::endl;
    id++;
  }
  if(Spectra::CompInfo::Successful != geigs.info()){
    std::cerr << "Failed " << __func__ <<std::endl;
    return false;
  }
    return true;
}

参考文献

1. Spectra, https://spectralib.org/
2. ブログページ「疎な巨大対称行列の一般化固有値問題のC++/Eigenでの効率的な解き方」　https://uncorrelated.hatenablog.com/entry/2020/10/02/060116