RokkoWrapperSerialDenseSolvers - t-sakashita/rokko GitHub Wiki

密行列ソルバ（逐次版）向けラッパーの実装

Eigen3の固有値ソルバ

http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html

  • lower partのみが使われる。
  • upper partを読むようには、できない。
  • Rokkoでは、eigen3でmatrix_part(またはuplow)がupperのとき、実行時エラーを出すようにした。

LAPACK

  • 現在、Rokkoでは、LAPACKのドライバルーチンを呼ぶようになっている。
  • Todo: 直接、三重対角化や三重対角行列向け固有値ルーチン(dstebzのsplitの対応など)を呼び出さないと使えない機能があるかも。

dsyevx

abstolの符号により、

  • 符号が正：２分法
  • 符号が負：QR法

ソースコード中のコメントより： http://www.netlib.org/lapack/explore-html/d2/d97/dsyevx_8f.html

If all eigenvalues are desired and ABSTOL is less than or equal to zero, then call DSTERF or DORGTR and SSTEQR. If this fails for some eigenvalue, then try DSTEBZ. Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.

Eigenvalues will be computed most accurately when ABSTOL is set to twice the underflow threshold 2DLAMCH('S'), not zero. If this routine returns with INFO>0, indicating that some eigenvectors did not converge, try setting ABSTOL to 2DLAMCH('S').

dsyevrの出力引数issupzの謎

This is relevant in the case when the matrix is split.

共通事項

以下を、固有値分解を行う行列、固有ベクトルを収める行列について別の値を設定することができる。

  • 固有ベクトル一本のサイズ
  • column-majorでは、行数。row-majorでは、列数。
  • LAPACKの場合、入力パラメータ'N'
  • Eigen3の場合
    • column-majorではrows()。row-majorではcols()
    • innerSize()で求められる。
  • 固有ベクトルの数
  • column-majorでは、列数。row-majorでは、行数。
  • LAPACKの場合、出力パラメータm
  • Eigen3の場合
    • column-majorではcols()。row-majorではrows()
    • outerSize()で求められる。

注意：上記と、stride（leading dimensionを含む）はまた別の概念。Eigen3のmajorやlocal dimensionについて

Todo

  • strideの次元と行列次元が一致しない場合（余白がある場合）も扱えるようにしたい。

Todo

  • verboseオプション

  • verboseレベルを用意

  • 表示する事項

  • 選択したソルバで使われなかったパラメータ

  • ソルバのヘルプ関数を用意

  • 必須パラメータ、オプションパラメータを出力

  • ソルバに渡す前に、入力パラメータを事前チェックする関数を用意。(LAPACKのエラー出力はわかりづらいため。)

  • LAPACK本体はFortranで書かれているため、行列はcol-majorである。LAPACKEをrow_majorで使うと、中で転置を行うことによるオーバーヘッドが起こる旨、警告する。

  • LAPACKに与える変数を、Rokkoラッパーでは構造体などでまとめた方が良いか。

  • LAPACKでは、lwork=-1によるサイズ問い合わせでのエラーはあり得るか？

備考

upper_valueとupper_index（lower_valueとlower_index）を共通にしなかった理由

当初、upper_limitとlower_limitで、以下の型によって自動的に意味を変えようとしていた。

  • int型：求める固有値の番号
  • double型：求める固有値の値の範囲
  • メリット：エラーチェックが楽になる。
  • デメリット：整数ではなく倍精度を意味したいのに、うっかりピリオドを忘れてしまうと、意味が変わってしまうため。
  params.set("upper_limit", 5.);  // upper_valueを意味
  params.set("upper_limit", 5);  // upper_indexを意味