概要
本サンプルはFortran言語によりLAPACKルーチンDGGESXを利用するサンプルプログラムです。
行列対

ここで


入力データ
(本ルーチンの詳細はDGGESX のマニュアルページを参照)1 2 3 4 5 6 7 8 9 10
このデータをダウンロード |
DGGESX 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
出力結果
(本ルーチンの詳細はDGGESX のマニュアルページを参照)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
この出力例をダウンロード |
DGGESX Example Program Results Matrix A 1 2 3 4 1 3.9000 12.5000 -34.5000 -0.5000 2 4.3000 21.5000 -47.5000 7.5000 3 4.3000 21.5000 -43.5000 3.5000 4 4.4000 26.0000 -46.0000 6.0000 Matrix B 1 2 3 4 1 1.0000 2.0000 -3.0000 1.0000 2 1.0000 3.0000 -5.0000 4.0000 3 1.0000 3.0000 -4.0000 3.0000 4 1.0000 3.0000 -4.0000 4.0000 Number of eigenvalues for which SELCTG is true = 2 (dimension of deflating subspaces) Selected generalized eigenvalues 1 ( 2.000, 0.000) 2 ( 4.000, 0.000) Reciprocals of left and right projection norms onto the deflating subspaces for the selected eigenvalues RCONDE(1) = 1.9E-01, RCONDE(2) = 1.8E-02 Reciprocal condition numbers for the left and right deflating subspaces RCONDV(1) = 5.4E-02, RCONDV(2) = 9.0E-02 Warning: Floating underflow occurred Approximate asymptotic error bound for selected eigenvalues = 1.1E-13 Approximate asymptotic error bound for the deflating subspaces = 2.4E-13
ソースコード
(本ルーチンの詳細はDGGESX のマニュアルページを参照)※本サンプルソースコードのご利用手順は「サンプルのコンパイル及び実行方法」をご参照下さい。
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 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
このソースコードをダウンロード |
! DGGESX Example Program Text ! Copyright 2017, Numerical Algorithms Group Ltd. http://www.nag.com Module dggesx_mod ! DGGESX Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use lapack_precision, Only: dp ! .. Implicit None Statement .. Implicit None ! .. Accessibility Statements .. Private Public :: selctg ! .. Parameters .. Integer, Parameter, Public :: nb = 64, nin = 5, nout = 6 Contains Function selctg(ar, ai, b) ! Logical function selctg for use with DGGESX (DGGESX) ! Returns the value .TRUE. if the eigenvalue is real and positive ! .. Function Return Value .. Logical :: selctg ! .. Scalar Arguments .. Real (Kind=dp), Intent (In) :: ai, ar, b ! .. Executable Statements .. selctg = (ar>0._dp .And. ai==0._dp .And. b/=0._dp) Return End Function End Module Program dggesx_example ! DGGESX Example Main Program ! .. Use Statements .. Use blas_interfaces, Only: dgemm Use dggesx_mod, Only: nb, nin, nout, selctg Use lapack_example_aux, Only: nagf_blas_dpyth, & nagf_file_print_matrix_real_gen Use lapack_interfaces, Only: dggesx, dlange Use lapack_precision, Only: dp ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=dp) :: abnorm, alph, anorm, bet, bnorm, eps, normd, norme, & tol Integer :: i, ifail, info, lda, ldb, ldc, ldd, lde, ldvsl, ldvsr, & liwork, lwork, n, sdim ! .. Local Arrays .. Real (Kind=dp), Allocatable :: a(:, :), alphai(:), alphar(:), b(:, :), & beta(:), c(:, :), d(:, :), e(:, :), vsl(:, :), vsr(:, :), work(:) Real (Kind=dp) :: rconde(2), rcondv(2), rdum(1) Integer :: idum(1) Integer, Allocatable :: iwork(:) Logical, Allocatable :: bwork(:) ! .. Intrinsic Procedures .. Intrinsic :: epsilon, max, nint ! .. Executable Statements .. Write (nout, *) 'DGGESX Example Program Results' Write (nout, *) Flush (nout) ! Skip heading in data file Read (nin, *) Read (nin, *) n lda = n ldb = n ldc = n ldd = n lde = n ldvsl = n ldvsr = n Allocate (a(lda,n), alphai(n), alphar(n), b(ldb,n), beta(n), & vsl(ldvsl,n), vsr(ldvsr,n), bwork(n), c(ldc,n), d(ldd,n), e(lde,n)) ! Use routine workspace query to get optimal workspace. lwork = -1 liwork = -1 Call dggesx('Vectors (left)', 'Vectors (right)', 'Sort', selctg, & 'Both reciprocal condition numbers', n, a, lda, b, ldb, sdim, alphar, & alphai, beta, vsl, ldvsl, vsr, ldvsr, rconde, rcondv, rdum, lwork, & idum, liwork, bwork, info) ! Make sure that there is enough workspace for block size nb. lwork = max(8*(n+1)+16+n*nb+n*n/2, nint(rdum(1))) liwork = max(n+6, idum(1)) Allocate (work(lwork), iwork(liwork)) ! 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) ! Copy A and B into D and E respectively d(1:n, 1:n) = a(1:n, 1:n) e(1:n, 1:n) = b(1:n, 1:n) ! Print matrices A and B ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call nagf_file_print_matrix_real_gen('General', ' ', n, n, a, lda, & 'Matrix A', ifail) Write (nout, *) Flush (nout) ifail = 0 Call nagf_file_print_matrix_real_gen('General', ' ', n, n, b, ldb, & 'Matrix B', ifail) Write (nout, *) Flush (nout) ! Find the Frobenius norms of A and B anorm = dlange('Frobenius', n, n, a, lda, work) bnorm = dlange('Frobenius', n, n, b, ldb, work) ! Find the generalized Schur form Call dggesx('Vectors (left)', 'Vectors (right)', 'Sort', selctg, & 'Both reciprocal condition numbers', n, a, lda, b, ldb, sdim, alphar, & alphai, beta, vsl, ldvsl, vsr, ldvsr, rconde, rcondv, work, lwork, & iwork, liwork, bwork, info) If (info==0 .Or. info==(n+2)) Then ! Compute A - Q*S*Z^T from the factorization of (A,B) and store in ! matrix D alph = 1.0_dp bet = 0.0_dp Call dgemm('N', 'N', n, n, n, alph, vsl, ldvsl, a, lda, bet, c, ldc) alph = -1.0_dp bet = 1.0_dp Call dgemm('N', 'T', n, n, n, alph, c, ldc, vsr, ldvsr, bet, d, ldd) ! Compute B - Q*T*Z^T from the factorization of (A,B) and store in ! matrix E alph = 1.0_dp bet = 0.0_dp Call dgemm('N', 'N', n, n, n, alph, vsl, ldvsl, b, ldb, bet, c, ldc) alph = -1.0_dp bet = 1.0_dp Call dgemm('N', 'T', n, n, n, alph, c, ldc, vsr, ldvsr, bet, e, lde) ! Find norms of matrices D and E and warn if either is too large normd = dlange('O', ldd, n, d, ldd, work) norme = dlange('O', lde, n, e, lde, work) If (normd>epsilon(1.0E0_dp)**0.8_dp .Or. norme>epsilon(1.0E0_dp)** & 0.8_dp) Then Write (nout, *) 'Norm of A-(Q*S*Z^T) or norm of B-(Q*T*Z^T) & &is much greater than 0.' Write (nout, *) 'Schur factorization has failed.' Else ! Print solution Write (nout, 100) & 'Number of eigenvalues for which SELCTG is true = ', sdim, & '(dimension of deflating subspaces)' Write (nout, *) ! Print generalized eigenvalues Write (nout, *) 'Selected generalized eigenvalues' Do i = 1, sdim If (beta(i)/=0.0_dp) Then Write (nout, 110) i, '(', alphar(i)/beta(i), ',', & alphai(i)/beta(i), ')' Else Write (nout, 120) i End If End Do If (info==(n+2)) Then Write (nout, 130) '***Note that rounding errors mean ', & 'that leading eigenvalues in the generalized', & 'Schur form no longer satisfy SELCTG = .TRUE.' Write (nout, *) End If Flush (nout) ! Print out the reciprocal condition numbers Write (nout, *) Write (nout, 140) & 'Reciprocals of left and right projection norms onto', & 'the deflating subspaces for the selected eigenvalues', & 'RCONDE(1) = ', rconde(1), ', RCONDE(2) = ', rconde(2) Write (nout, *) Write (nout, 140) & 'Reciprocal condition numbers for the left and right', & 'deflating subspaces', 'RCONDV(1) = ', rcondv(1), & ', RCONDV(2) = ', rcondv(2) Flush (nout) ! Compute the machine precision and sqrt(anorm**2+bnorm**2) eps = epsilon(1.0E0_dp) abnorm = nagf_blas_dpyth(anorm, bnorm) tol = eps*abnorm ! Print out the approximate asymptotic error bound on the ! average absolute error of the selected eigenvalues given by ! eps*norm((A, B))/PL, where PL = RCONDE(1) Write (nout, *) Write (nout, 150) 'Approximate asymptotic error bound for selected ' & , 'eigenvalues = ', tol/rconde(1) ! Print out an approximate asymptotic bound on the maximum ! angular error in the computed deflating subspaces given by ! eps*norm((A, B))/DIF(2), where DIF(2) = RCONDV(2) Write (nout, 150) & 'Approximate asymptotic error bound for the deflating ', & 'subspaces = ', tol/rcondv(2) End If Else Write (nout, 100) 'Failure in DGGESX. INFO =', info End If 100 Format (1X, A, I4, /, 1X, A) 110 Format (1X, I4, 5X, A, F7.3, A, F7.3, A) 120 Format (1X, I4, 'Eigenvalue is infinite') 130 Format (1X, 2A, /, 1X, A) 140 Format (1X, A, /, 1X, A, /, 1X, 2(A,1P,E8.1)) 150 Format (1X, 2A, 1P, E8.1) End Program