fortran90でLapackを用います。

倍精度ルーチン

・複素一般行列の対角化, diag(N,A,Ev)
・複素エルミート行列の対角化, zdiag(N,A,Ev)
・実エルミート行列の対角化, herdiag(N,A,Ev), diag_slct(N,A,Nq,Ev)
・実対称三重対角行列の対角化, tridiag(N,A,Ev)

4倍精度ルーチン

・複素一般行列の対角化, qdiag(N,A,Ev)

複素一般行列　　　　　　　　　：行列要素は複素数で、値は適当に入っている。
複素エルミート行列　　　　　　：行列要素は複素数で、\(A_{i,j}=A^*_{j,i}\)の関係がある。
実エルミート行列（実対称行列）：行列要素は実数で、\(A_{i,j}=A_{j,i}\)の関係がある。
実三重対角行列　　　　　　　　：行列要素は実数で、\(A_{i,i}\ne 0, ~~A_{i,i-1}=A_{i,i+1}\)、それ以外の要素は0。

（エルミートな行列と言われたら、その固有値は必ず実数です。）

複素一般行列の場合

wolframのページで、
Eigenvectors[{-2,2,4},{-2,4,2},{-2,1,4}]
と入力することで固有値、固有ベクトルを出してくれます。これで確かめを行いました。

・この行列を入れた時、正しい固有値と固有ベクトルを出力してくれる事を確認しています。
・40×40程度の行列でも確認しています。
・複素数の固有値、固有ベクトルでも確認しています。

Lapackを使います。

か、g95, もしくはgfortranでlapackにリンク

して、コンパイルしてください。

大きさN×Nの正方行列Aを入力すると, 固有値Ev(1:N)とそれぞれの固有値の要素に属する固有ベクトルを上書きして返します。(1列目に固有値Ev(1)に属する固有ベクトルが、2列目に固有値Ev(2)に属する固有ベクトルが･･･という感じに。)
また固有値は実部の大小を比較して小さい順に並べ替えられ、固有ベクトルもそれに対応した列に格納されます。
計算速度は\(O(N^3)\)です。

複素エルミート行列の対角化

対称行列に限る場合は計算を早くすることができます。
ここでは複素エルミート行列を対角化して固有値を出力するケースを考えましょう。

そのまま使うのはやはり面倒です。なぜなら、このルーチンを用いるためにはワーク配列を何個か宣言して用いる必要があるからです。
このワーク配列は計算に直接関係ありません。
僕らが欲しいのは行列を入れたら固有値と固有ベクトルを返してくれさえすればいいものです。
下はサブルーチンzdiagは僕が作ったもので、複素エルミート行列\(A(1:N,1:N)\)を入れると固有値の配列\(Ev(1:N)\)と\(A\)に上書きして固有ベクトル\(A(1:N,1:N)\)を入れて返します。

エルミート行列なので、実は全ての行列要素を求める必要ありません。
下の例では、
エルミート行列\(A(i,j)\)は,\(j\ge i\)の上三角の要素だけ分かっていればいいです。
なので、コメントアウトしている行列でも同じ固有値と固有ベクトルが得られます。
以下のプログラムをlapackと共にコンパイルし、実行すると、

という結果が得られます。
1番目の固有値(-0.42419E+01)に対応する固有ベクトルは1列目に、
2番目の固有値( 0.21580E+01)に対応する固有ベクトルは2列目…というように代入されていきます。

プログラム例)
確かめ用のリンクEigenvectors[{-2,2+i,4},{2-i,4,1-2i},{4,1+2i,4}] -wolfram alpha

行と列を入れ替えたければ、以下のルーチンを使用してください。

実エルミート行列の場合

以下のサブルーチンを使うと良いでしょう。

\(N\times N\)の実対称行列で、固有値の小さい方から順に\(Nq\)個だけ欲しい場合

実エルミート行列の、
完全対角化(herdiag)と、
選択的に固有値を計算する(diag_slct)
の比較(実時間計測)を行いました。

横軸が選択的に何個固有値を計算するか、を表し、縦軸は実際に掛かった時間をあらわします。
色は行列サイズによって変化し、黄色,紫、緑の順に行列サイズが\(2^{16}\times 2^{16}, 2^{14}\times 2^{14}, 2^{12}\times 2^{12}\)を表します。
破線はherdiagによる対角化時間、実線はdiag_slctによる計算時間を表します。
また、図中の値2,4,8は計算に用いたCPUの数です。CPUは8コア16スレッドです。

ここから分かるのは、欲しい固有値の数が行列サイズ\(N\times N\)の\(N/4\)より小さい場合は選択的に計算した方が早い、ということです。
選択的に計算する方は計算時に要求されるメモリも少なくて済みます。

(2016/08/27)

4倍精度の一般行列の対角化

4倍精度用の一般行列(複素非エルミートを含む)を対角化し、その固有値と固有ベクトルを出力するfortranコードが以下のページ
４倍精度対角化 -ただの備忘録
で公開されていました。
ここで公開されているサブルーチン”qeispack.f90”
を使うことで対角化できます。
gfortranとifortで動くことは確かめました。
コンパイルは例えば

のようにしてください。
本ページと同じように書き加えて実装すると以下のようになります。

実対称三重対角行列の場合

実対称三角行列の場合も置いておきます。

Tweet