LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zget52.f
Go to the documentation of this file.
1 *> \brief \b ZGET52
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
12 * WORK, RWORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL LEFT
16 * INTEGER LDA, LDB, LDE, N
17 * ..
18 * .. Array Arguments ..
19 * DOUBLE PRECISION RESULT( 2 ), RWORK( * )
20 * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
21 * $ BETA( * ), E( LDE, * ), WORK( * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> ZGET52 does an eigenvector check for the generalized eigenvalue
31 *> problem.
32 *>
33 *> The basic test for right eigenvectors is:
34 *>
35 *> | b(i) A E(i) - a(i) B E(i) |
36 *> RESULT(1) = max -------------------------------
37 *> i n ulp max( |b(i) A|, |a(i) B| )
38 *>
39 *> using the 1-norm. Here, a(i)/b(i) = w is the i-th generalized
40 *> eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
41 *> generalized eigenvalue of m A - B.
42 *>
43 *> H H _ _
44 *> For left eigenvectors, A , B , a, and b are used.
45 *>
46 *> ZGET52 also tests the normalization of E. Each eigenvector is
47 *> supposed to be normalized so that the maximum "absolute value"
48 *> of its elements is 1, where in this case, "absolute value"
49 *> of a complex value x is |Re(x)| + |Im(x)| ; let us call this
50 *> maximum "absolute value" norm of a vector v M(v).
51 *> If a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
52 *> vector. The normalization test is:
53 *>
54 *> RESULT(2) = max | M(v(i)) - 1 | / ( n ulp )
55 *> eigenvectors v(i)
56 *>
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] LEFT
63 *> \verbatim
64 *> LEFT is LOGICAL
65 *> =.TRUE.: The eigenvectors in the columns of E are assumed
66 *> to be *left* eigenvectors.
67 *> =.FALSE.: The eigenvectors in the columns of E are assumed
68 *> to be *right* eigenvectors.
69 *> \endverbatim
70 *>
71 *> \param[in] N
72 *> \verbatim
73 *> N is INTEGER
74 *> The size of the matrices. If it is zero, ZGET52 does
75 *> nothing. It must be at least zero.
76 *> \endverbatim
77 *>
78 *> \param[in] A
79 *> \verbatim
80 *> A is COMPLEX*16 array, dimension (LDA, N)
81 *> The matrix A.
82 *> \endverbatim
83 *>
84 *> \param[in] LDA
85 *> \verbatim
86 *> LDA is INTEGER
87 *> The leading dimension of A. It must be at least 1
88 *> and at least N.
89 *> \endverbatim
90 *>
91 *> \param[in] B
92 *> \verbatim
93 *> B is COMPLEX*16 array, dimension (LDB, N)
94 *> The matrix B.
95 *> \endverbatim
96 *>
97 *> \param[in] LDB
98 *> \verbatim
99 *> LDB is INTEGER
100 *> The leading dimension of B. It must be at least 1
101 *> and at least N.
102 *> \endverbatim
103 *>
104 *> \param[in] E
105 *> \verbatim
106 *> E is COMPLEX*16 array, dimension (LDE, N)
107 *> The matrix of eigenvectors. It must be O( 1 ).
108 *> \endverbatim
109 *>
110 *> \param[in] LDE
111 *> \verbatim
112 *> LDE is INTEGER
113 *> The leading dimension of E. It must be at least 1 and at
114 *> least N.
115 *> \endverbatim
116 *>
117 *> \param[in] ALPHA
118 *> \verbatim
119 *> ALPHA is COMPLEX*16 array, dimension (N)
120 *> The values a(i) as described above, which, along with b(i),
121 *> define the generalized eigenvalues.
122 *> \endverbatim
123 *>
124 *> \param[in] BETA
125 *> \verbatim
126 *> BETA is COMPLEX*16 array, dimension (N)
127 *> The values b(i) as described above, which, along with a(i),
128 *> define the generalized eigenvalues.
129 *> \endverbatim
130 *>
131 *> \param[out] WORK
132 *> \verbatim
133 *> WORK is COMPLEX*16 array, dimension (N**2)
134 *> \endverbatim
135 *>
136 *> \param[out] RWORK
137 *> \verbatim
138 *> RWORK is DOUBLE PRECISION array, dimension (N)
139 *> \endverbatim
140 *>
141 *> \param[out] RESULT
142 *> \verbatim
143 *> RESULT is DOUBLE PRECISION array, dimension (2)
144 *> The values computed by the test described above. If A E or
145 *> B E is likely to overflow, then RESULT(1:2) is set to
146 *> 10 / ulp.
147 *> \endverbatim
148 *
149 * Authors:
150 * ========
151 *
152 *> \author Univ. of Tennessee
153 *> \author Univ. of California Berkeley
154 *> \author Univ. of Colorado Denver
155 *> \author NAG Ltd.
156 *
157 *> \date November 2011
158 *
159 *> \ingroup complex16_eig
160 *
161 * =====================================================================
162  SUBROUTINE zget52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
163  $ work, rwork, result )
164 *
165 * -- LAPACK test routine (version 3.4.0) --
166 * -- LAPACK is a software package provided by Univ. of Tennessee, --
167 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168 * November 2011
169 *
170 * .. Scalar Arguments ..
171  LOGICAL left
172  INTEGER lda, ldb, lde, n
173 * ..
174 * .. Array Arguments ..
175  DOUBLE PRECISION result( 2 ), rwork( * )
176  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
177  $ beta( * ), e( lde, * ), work( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  DOUBLE PRECISION zero, one
184  parameter( zero = 0.0d+0, one = 1.0d+0 )
185  COMPLEX*16 czero, cone
186  parameter( czero = ( 0.0d+0, 0.0d+0 ),
187  $ cone = ( 1.0d+0, 0.0d+0 ) )
188 * ..
189 * .. Local Scalars ..
190  CHARACTER normab, trans
191  INTEGER j, jvec
192  DOUBLE PRECISION abmax, alfmax, anorm, betmax, bnorm, enorm,
193  $ enrmer, errnrm, safmax, safmin, scale, temp1,
194  $ ulp
195  COMPLEX*16 acoeff, alphai, bcoeff, betai, x
196 * ..
197 * .. External Functions ..
198  DOUBLE PRECISION dlamch, zlange
199  EXTERNAL dlamch, zlange
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL zgemv
203 * ..
204 * .. Intrinsic Functions ..
205  INTRINSIC abs, dble, dconjg, dimag, max
206 * ..
207 * .. Statement Functions ..
208  DOUBLE PRECISION abs1
209 * ..
210 * .. Statement Function definitions ..
211  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
212 * ..
213 * .. Executable Statements ..
214 *
215  result( 1 ) = zero
216  result( 2 ) = zero
217  IF( n.LE.0 )
218  $ return
219 *
220  safmin = dlamch( 'Safe minimum' )
221  safmax = one / safmin
222  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
223 *
224  IF( left ) THEN
225  trans = 'C'
226  normab = 'I'
227  ELSE
228  trans = 'N'
229  normab = 'O'
230  END IF
231 *
232 * Norm of A, B, and E:
233 *
234  anorm = max( zlange( normab, n, n, a, lda, rwork ), safmin )
235  bnorm = max( zlange( normab, n, n, b, ldb, rwork ), safmin )
236  enorm = max( zlange( 'O', n, n, e, lde, rwork ), ulp )
237  alfmax = safmax / max( one, bnorm )
238  betmax = safmax / max( one, anorm )
239 *
240 * Compute error matrix.
241 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
242 *
243  DO 10 jvec = 1, n
244  alphai = alpha( jvec )
245  betai = beta( jvec )
246  abmax = max( abs1( alphai ), abs1( betai ) )
247  IF( abs1( alphai ).GT.alfmax .OR. abs1( betai ).GT.betmax .OR.
248  $ abmax.LT.one ) THEN
249  scale = one / max( abmax, safmin )
250  alphai = scale*alphai
251  betai = scale*betai
252  END IF
253  scale = one / max( abs1( alphai )*bnorm, abs1( betai )*anorm,
254  $ safmin )
255  acoeff = scale*betai
256  bcoeff = scale*alphai
257  IF( left ) THEN
258  acoeff = dconjg( acoeff )
259  bcoeff = dconjg( bcoeff )
260  END IF
261  CALL zgemv( trans, n, n, acoeff, a, lda, e( 1, jvec ), 1,
262  $ czero, work( n*( jvec-1 )+1 ), 1 )
263  CALL zgemv( trans, n, n, -bcoeff, b, lda, e( 1, jvec ), 1,
264  $ cone, work( n*( jvec-1 )+1 ), 1 )
265  10 continue
266 *
267  errnrm = zlange( 'One', n, n, work, n, rwork ) / enorm
268 *
269 * Compute RESULT(1)
270 *
271  result( 1 ) = errnrm / ulp
272 *
273 * Normalization of E:
274 *
275  enrmer = zero
276  DO 30 jvec = 1, n
277  temp1 = zero
278  DO 20 j = 1, n
279  temp1 = max( temp1, abs1( e( j, jvec ) ) )
280  20 continue
281  enrmer = max( enrmer, temp1-one )
282  30 continue
283 *
284 * Compute RESULT(2) : the normalization error in E.
285 *
286  result( 2 ) = enrmer / ( dble( n )*ulp )
287 *
288  return
289 *
290 * End of ZGET52
291 *
292  END