概要
本サンプルはFortran言語によりLAPACKルーチンDGGEVXを利用するサンプルプログラムです。
行列対
の全ての固有値と右固有ベクトルを求めます。
条件数の推定値とそれぞれの固有値と固有ベクトルの前方誤差限界も合わせて求めます。行列対を均衡化するオプションが指定されます。
入力データ
（本ルーチンの詳細はDGGEVX のマニュアルページを参照）1 2 3 4 5 6 7 8 9 10
| このデータをダウンロード |
DGGEVX 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
出力結果
（本ルーチンの詳細はDGGEVX のマニュアルページを参照）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
| この出力例をダウンロード |
Warning: Floating underflow occurred
DGGEVX Example Program Results
Eigenvalue( 1) = 2.00000
Eigenvector( 1)
0.99606
0.00569
0.06261
0.06261
Eigenvalue( 2) = ( 3.00000, -4.00000)
Eigenvector( 2)
( 0.94491, -0.00000)
( 0.18898, 0.00000)
( 0.11339, 0.15119)
( 0.11339, 0.15119)
Eigenvalue( 3) = ( 3.00000, 4.00000)
Eigenvector( 3)
( 0.94491, 0.00000)
( 0.18898, -0.00000)
( 0.11339, -0.15119)
( 0.11339, -0.15119)
Eigenvalue( 4) = 4.00000
Eigenvector( 4)
0.98752
0.01097
-0.03292
0.15361
ソースコード
（本ルーチンの詳細はDGGEVX のマニュアルページを参照）※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
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 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
| このソースコードをダウンロード |
Program dggevx_example
! DGGEVX Example Program Text
! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com
! .. Use Statements ..
Use blas_interfaces, Only: dnrm2
Use lapack_example_aux, Only: nagf_sort_realvec_rank, &
nagf_sort_realvec_rank_rearrange
Use lapack_interfaces, Only: dggevx
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nb = 64, nin = 5, nout = 6
Logical, Parameter :: verbose = .False.
! .. Local Scalars ..
Complex (Kind=dp) :: eig
Real (Kind=dp) :: abnrm, bbnrm, eps, jswap, rcnd, scal_i, scal_r, small
Integer :: i, ifail, ihi, ilo, info, j, k, lda, ldb, ldvr, lwork, n
Logical :: pair
! .. Local Arrays ..
Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), &
beta(:), lscale(:), rconde(:), rcondv(:), rscale(:), vr(:, :), &
vr_row(:), work(:)
Real (Kind=dp) :: dummy(1, 1)
Integer, Allocatable :: iwork(:)
Logical, Allocatable :: bwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, all, cmplx, epsilon, max, maxloc, nint, sqrt, sum, &
tiny
! .. Executable Statements ..
Write (nout, *) 'DGGEVX 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), lscale(n), &
rconde(n), rcondv(n), rscale(n), vr(ldvr,n), iwork(n+6), bwork(n))
! Use routine workspace query to get optimal workspace.
lwork = -1
Call dggevx('Balance', 'No vectors (left)', 'Vectors (right)', &
'Both reciprocal condition numbers', n, a, lda, b, ldb, alphar, &
alphai, beta, dummy, 1, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, &
bbnrm, rconde, rcondv, dummy, lwork, iwork, bwork, info)
! Make sure that there is enough workspace for block size nb.
lwork = max((nb+2*n)*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 dggevx('Balance', 'No vectors (left)', 'Vectors (right)', &
'Both reciprocal condition numbers', n, a, lda, b, ldb, alphar, &
alphai, beta, dummy, 1, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, &
bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
If (info>0) Then
Write (nout, *)
Write (nout, 100) 'Failure in DGGEVX. INFO =', info
Else
! Compute the machine precision and the safe range parameter
! small
eps = epsilon(1.0E0_dp)
small = tiny(1.0E0_dp)
! If beta(:) > eps, Order eigenvalues by ascending real parts
If (all(abs(beta(1:n))>eps)) Then
work(1:n) = alphar(1:n)/beta(1:n)
ifail = 0
Call nagf_sort_realvec_rank(work, 1, n, 'Ascending', iwork, ifail)
Call nagf_sort_realvec_rank_rearrange(alphar, 1, n, iwork, ifail)
Call nagf_sort_realvec_rank_rearrange(alphai, 1, n, iwork, ifail)
Call nagf_sort_realvec_rank_rearrange(beta, 1, n, iwork, 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, iwork, ifail)
vr(j, 1:n) = vr_row(1:n)
End Do
Deallocate (vr_row)
End If
! Print out eigenvalues and vectors and associated condition
! number and bounds
pair = .False.
Do j = 1, n
! Print out information on the j-th eigenvalue
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 (.Not. pair) Then
jswap = 1.0_dp
If (alphai(j)>0.0_dp) Then
jswap = -jswap
End If
End If
If (alphai(j)==0.0E0_dp) Then
Write (nout, 120) 'Eigenvalue(', j, ') = ', alphar(j)/beta(j)
Else
eig = cmplx(alphar(j), jswap*alphai(j), kind=dp)/ &
cmplx(beta(j), kind=dp)
Write (nout, 130) 'Eigenvalue(', j, ') = ', eig
End If
End If
If (verbose) Then
rcnd = rconde(j)
Write (nout, *)
Write (nout, 140) ' Reciprocal condition number = ', rcnd
End If
! Print out information on the j-th eigenvector
Write (nout, *)
Write (nout, 150) 'Eigenvector(', j, ')'
If (alphai(j)==0.0E0_dp) Then
! Let largest element be positive
work(1:n) = abs(vr(1:n,j))
k = maxloc(work(1:n), 1)
If (vr(k,j)<0.0_dp) Then
vr(1:n, j) = -vr(1:n, j)/dnrm2(n, vr(1,j), 1)
End If
Write (nout, 160)(vr(i,j), i=1, n)
Else
If (pair) Then
Write (nout, 170)(vr(i,j-1), -jswap*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))/sqrt(sum(work(1:n)))
scal_i = -vr(k, j+1)/sqrt(work(k))/sqrt(sum(work(1:n)))
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)
vr(k, j+1) = 0.0_dp
Write (nout, 170)(vr(i,j), jswap*vr(i,j+1), i=1, n)
End If
pair = .Not. pair
End If
If (verbose) Then
rcnd = rcondv(j)
Write (nout, *)
Write (nout, 140) ' Reciprocal condition number = ', rcnd
End If
End Do
End If
100 Format (1X, A, I4)
110 Format (/, 1X, A, I2, 2A, /, 1X, 2(A,I2,A,F11.5), A, I2, A, F11.5)
120 Format (/, 1X, A, I2, A, F11.5)
130 Format (/, 1X, A, I2, A, '(', F11.5, ',', F11.5, ')')
140 Format (1X, A, 1P, E8.1)
150 Format (1X, A, I2, A)
160 Format (1X, F11.5)
170 Format (1X, '(', F11.5, ',', F11.5, ')')
End Program