概要
本サンプルはFortran言語によりLAPACKルーチンDGGEVを利用するサンプルプログラムです。
行列対
の全ての固有値と右固有ベクトルを求めます。
入力データ
（本ルーチンの詳細はDGGEV のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10
| このデータをダウンロード |
DGGEV Example Program Data 4 :Value of N 3.9 12.5 -34.5 -0.5 4.3 21.5 -47.5 7.5 4.3 21.5 -43.5 3.5 4.4 26.0 -46.0 6.0 :End of matrix A 1.0 2.0 -3.0 1.0 1.0 3.0 -5.0 4.0 1.0 3.0 -4.0 3.0 1.0 3.0 -4.0 4.0 :End of matrix B
出力結果
（本ルーチンの詳細はDGGEV のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
| この出力例をダウンロード |
Warning: Floating underflow occurred
DGGEV Example Program Results
Eigenvalue( 1) = 20.000
Eigenvector( 1)
10.00000
0.05714
0.62857
0.62857
Eigenvalue( 2) = ( 30.000, -40.000)
Eigenvector( 2)
( 7.12215, 0.00000)
( 1.42443, -0.00000)
( 0.85466, 1.13954)
( 0.85466, 1.13954)
Eigenvalue( 3) = ( 30.000, 40.000)
Eigenvector( 3)
( 7.12215, -0.00000)
( 1.42443, 0.00000)
( 0.85466, -1.13954)
( 0.85466, -1.13954)
Eigenvalue( 4) = 40.000
Eigenvector( 4)
10.00000
0.11111
-0.33333
1.55556
ソースコード
（本ルーチンの詳細はDGGEV のマニュアルページを参照）※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
| このソースコードをダウンロード |
Program dggev_example
! DGGEV Example Program Text
! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com
! .. Use Statements ..
Use lapack_example_aux, Only: nagf_sort_realmat_rank_rows, &
nagf_sort_realvec_rank_rearrange
Use lapack_interfaces, Only: dggev
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) :: scal_i, scal_r, small
Integer :: i, ifail, info, j, k, lda, ldb, ldvr, lwork, n
Logical :: pair
! .. Local Arrays ..
Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
beta(:), vr(:, :), vr_row(:), work(:)
Real (Kind=dp) :: dummy(1, 1)
Integer, Allocatable :: irank(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, all, cmplx, epsilon, max, maxloc, nint, sqrt, tiny
! .. Executable Statements ..
Write (nout, *) 'DGGEV Example Program Results'
! Skip heading in data file
Read (nin, *)
Read (nin, *) n
lda = n
ldb = n
ldvr = n
Allocate (a(lda,n), alphai(n), alphar(n), b(ldb,n), beta(n), vr(ldvr,n), &
irank(n))
! Use routine workspace query to get optimal workspace.
lwork = -1
Call dggev('No left vectors', 'Vectors (right)', n, a, lda, b, ldb, &
alphar, alphai, beta, dummy, 1, vr, ldvr, dummy, lwork, info)
! Make sure that there is enough workspace for block size nb.
lwork = max((nb+7)*n, nint(dummy(1,1)))
Allocate (work(lwork))
! Read in the matrices A and B
Read (nin, *)(a(i,1:n), i=1, n)
Read (nin, *)(b(i,1:n), i=1, n)
! Solve the generalized eigenvalue problem
Call dggev('No left vectors', 'Vectors (right)', n, a, lda, b, ldb, &
alphar, alphai, beta, dummy, 1, vr, ldvr, work, lwork, info)
If (info>0) Then
Write (nout, *)
Write (nout, 100) 'Failure in DGGEV. INFO =', info
Else
! If beta(:) > eps, Order eigenvalues by ascending real parts
! and then by ascending imaginary parts
If (all(abs(beta(1:n))>epsilon(1.0E0_dp))) Then
work(1:n) = alphar(1:n)/beta(1:n)
work(n+1:2*n) = alphai(1:n)/beta(1:n)
ifail = 0
Call nagf_sort_realmat_rank_rows(work, n, 1, n, 1, 2, 'Ascending', &
irank, ifail)
Call nagf_sort_realvec_rank_rearrange(alphar, 1, n, irank, ifail)
Call nagf_sort_realvec_rank_rearrange(alphai, 1, n, irank, ifail)
Call nagf_sort_realvec_rank_rearrange(beta, 1, n, irank, ifail)
! Order the eigenvectors in the same way
Allocate (vr_row(n))
Do j = 1, n
vr_row(1:n) = vr(j, 1:n)
Call nagf_sort_realvec_rank_rearrange(vr_row, 1, n, irank, ifail)
vr(j, 1:n) = vr_row(1:n)
End Do
Deallocate (vr_row)
End If
small = tiny(1.0E0_dp)
pair = .False.
Do j = 1, n
Write (nout, *)
If ((abs(alphar(j))+abs(alphai(j)))*small>=abs(beta(j))) Then
Write (nout, 110) 'Eigenvalue(', j, ')', &
' is numerically infinite or undetermined', 'ALPHAR(', j, &
') = ', alphar(j), ', ALPHAI(', j, ') = ', alphai(j), ', BETA(', &
j, ') = ', beta(j)
Else
If (alphai(j)==zero) Then
Write (nout, 120) 'Eigenvalue(', j, ') = ', alphar(j)/beta(j)
Else
eig = cmplx(alphar(j), alphai(j), kind=dp)/ &
cmplx(beta(j), kind=dp)
Write (nout, 130) 'Eigenvalue(', j, ') = ', eig
End If
End If
Write (nout, *)
Write (nout, 140) 'Eigenvector(', j, ')'
If (alphai(j)==zero) Then
! Let largest element be positive
work(1:n) = abs(vr(1:n,j))
k = maxloc(work(1:n), 1)
If (vr(k,j)<zero) Then
vr(1:n, j) = -vr(1:n, j)
End If
Write (nout, 150)(vr(i,j), i=1, n)
Else
If (pair) Then
Write (nout, 160)(vr(i,j-1), -vr(i,j), i=1, n)
Else
! Let largest element be real (and positive).
work(1:n) = vr(1:n, j)**2 + vr(1:n, j+1)**2
k = maxloc(work(1:n), 1)
scal_r = vr(k, j)/sqrt(work(k))
scal_i = -vr(k, j+1)/sqrt(work(k))
work(1:n) = vr(1:n, j)
vr(1:n, j) = scal_r*work(1:n) - scal_i*vr(1:n, j+1)
vr(1:n, j+1) = scal_r*vr(1:n, j+1) + scal_i*work(1:n)
Write (nout, 160)(vr(i,j), vr(i,j+1), i=1, n)
End If
pair = .Not. pair
End If
End Do
End If
100 Format (1X, A, I4)
110 Format (1X, A, I2, 2A, /, 1X, 2(A,I2,A,1P,F11.3,3X), A, I2, A, 1P, &
F11.3)
120 Format (1X, A, I2, A, 1P, F11.3)
130 Format (1X, A, I2, A, '(', 1P, F11.3, ',', 1P, F11.3, ')')
140 Format (1X, A, I2, A)
150 Format (1X, 1P, F11.5)
160 Format (1X, '(', 1P, F11.5, ',', 1P, F11.5, ')')
End Program