概要
本サンプルはFortran言語によりLAPACKルーチンDGEEVを利用するサンプルプログラムです。
行列のすべての固有値と右固有ベクトルを求めます。
入力データ
(本ルーチンの詳細はDGEEV のマニュアルページを参照)このデータをダウンロード |
DGEEV Example Program Data 4 :Value of N 0.35 0.45 -0.14 -0.17 0.09 0.07 -0.54 0.35 -0.44 -0.33 -0.03 0.17 0.25 -0.32 -0.13 0.11 :End of matrix A
出力結果
(本ルーチンの詳細はDGEEV のマニュアルページを参照)この出力例をダウンロード |
DGEEV Example Program Results Eigenvalue( 1) = 7.9948E-01 Eigenvector( 1) -6.5509E-01 -5.2363E-01 5.3622E-01 -9.5607E-02 Eigenvalue( 2) = (-9.9412E-02, 4.0079E-01) Eigenvector( 2) (-1.9330E-01, 2.5463E-01) ( 2.5186E-01,-5.2240E-01) ( 9.7182E-02,-3.0838E-01) ( 6.7595E-01, 0.0000E+00) Eigenvalue( 3) = (-9.9412E-02,-4.0079E-01) Eigenvector( 3) (-1.9330E-01,-2.5463E-01) ( 2.5186E-01, 5.2240E-01) ( 9.7182E-02, 3.0838E-01) ( 6.7595E-01,-0.0000E+00) Eigenvalue( 4) = -1.0066E-01 Eigenvector( 4) 1.2533E-01 3.3202E-01 5.9384E-01 7.2209E-01
ソースコード
(本ルーチンの詳細はDGEEV のマニュアルページを参照)※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
このソースコードをダウンロード |
Program dgeev_example ! DGEEV Example Program Text ! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com ! .. Use Statements .. Use lapack_interfaces, Only: dgeev Use lapack_precision, Only: dp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=dp), Parameter :: zero = 0.0_dp Integer, Parameter :: nb = 64, nin = 5, nout = 6 ! .. Local Scalars .. Complex (Kind=dp) :: eig Real (Kind=dp) :: alpha, beta, scale Integer :: i, info, j, k, lda, ldvr, lwork, n ! .. Local Arrays .. Real (Kind=dp), Allocatable :: a(:, :), vr(:, :), wi(:), work(:), wr(:) Real (Kind=dp) :: dummy(1, 1) ! .. Intrinsic Procedures .. Intrinsic :: cmplx, max, maxloc, nint, sqrt ! .. Executable Statements .. Write (nout, *) 'DGEEV Example Program Results' ! Skip heading in data file Read (nin, *) Read (nin, *) n lda = n ldvr = n Allocate (a(lda,n), vr(ldvr,n), wi(n), wr(n)) ! Use routine workspace query to get optimal workspace. lwork = -1 Call dgeev('No left vectors', 'Vectors (right)', n, a, lda, wr, wi, & dummy, 1, vr, ldvr, dummy, lwork, info) ! Make sure that there is enough workspace for block size nb. lwork = max((nb+2)*n, nint(dummy(1,1))) Allocate (work(lwork)) ! Read the matrix A from data file Read (nin, *)(a(i,1:n), i=1, n) ! Compute the eigenvalues and right eigenvectors of A Call dgeev('No left vectors', 'Vectors (right)', n, a, lda, wr, wi, & dummy, 1, vr, ldvr, work, lwork, info) If (info==0) Then ! Print solution Do j = 1, n Write (nout, *) If (wi(j)==zero) Then Write (nout, 100) 'Eigenvalue(', j, ') = ', wr(j) Else eig = cmplx(wr(j), wi(j), kind=dp) Write (nout, 110) 'Eigenvalue(', j, ') = ', eig End If Write (nout, *) Write (nout, 120) 'Eigenvector(', j, ')' If (wi(j)==zero) Then ! Scale by making largest element positive k = maxloc(vr(1:n,j), 1) If (vr(k,j)<zero) Then vr(1:n, j) = -vr(1:n, j) End If Write (nout, 130) vr(1:n, j) Else If (wi(j)>0.0E0_dp) Then ! Scale by making largest element real and positive work(1:n) = vr(1:n, j)**2 + vr(1:n, j+1)**2 k = maxloc(work(1:n), 1) scale = sqrt(work(k)) work(1:n) = vr(1:n, j) alpha = vr(k, j)/scale beta = vr(k, j+1)/scale vr(1:n, j) = alpha*work(1:n) + beta*vr(1:n, j+1) vr(1:n, j+1) = alpha*vr(1:n, j+1) - beta*work(1:n) Write (nout, 140)(vr(i,j), vr(i,j+1), i=1, n) Else Write (nout, 140)(vr(i,j-1), -vr(i,j), i=1, n) End If End Do Else Write (nout, *) Write (nout, 150) 'Failure in DGEEV. INFO = ', info End If 100 Format (1X, A, I2, A, 1P, E11.4) 110 Format (1X, A, I2, A, '(', 1P, E11.4, ',', 1P, E11.4, ')') 120 Format (1X, A, I2, A) 130 Format (1X, 1P, E11.4) 140 Format (1X, '(', 1P, E11.4, ',', 1P, E11.4, ')') 150 Format (1X, A, I4) End Program