LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dget52.f
Go to the documentation of this file.
1 *> \brief \b DGET52
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 DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
12 * ALPHAI, BETA, WORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL LEFT
16 * INTEGER LDA, LDB, LDE, N
17 * ..
18 * .. Array Arguments ..
19 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
20 * $ B( LDB, * ), BETA( * ), E( LDE, * ),
21 * $ RESULT( 2 ), WORK( * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> DGET52 does an eigenvector check for the generalized eigenvalue
31 *> problem.
32 *>
33 *> The basic test for right eigenvectors is:
34 *>
35 *> | b(j) A E(j) - a(j) B E(j) |
36 *> RESULT(1) = max -------------------------------
37 *> j n ulp max( |b(j) A|, |a(j) B| )
38 *>
39 *> using the 1-norm. Here, a(j)/b(j) = w is the j-th generalized
40 *> eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
41 *> generalized eigenvalue of m A - B.
42 *>
43 *> For real eigenvalues, the test is straightforward. For complex
44 *> eigenvalues, E(j) and a(j) are complex, represented by
45 *> Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
46 *> eigenvector becomes
47 *>
48 *> max( |Wr|, |Wi| )
49 *> --------------------------------------------
50 *> n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
51 *>
52 *> where
53 *>
54 *> Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
55 *>
56 *> Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
57 *>
58 *> T T _
59 *> For left eigenvectors, A , B , a, and b are used.
60 *>
61 *> DGET52 also tests the normalization of E. Each eigenvector is
62 *> supposed to be normalized so that the maximum "absolute value"
63 *> of its elements is 1, where in this case, "absolute value"
64 *> of a complex value x is |Re(x)| + |Im(x)| ; let us call this
65 *> maximum "absolute value" norm of a vector v M(v).
66 *> if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
67 *> vector. The normalization test is:
68 *>
69 *> RESULT(2) = max | M(v(j)) - 1 | / ( n ulp )
70 *> eigenvectors v(j)
71 *> \endverbatim
72 *
73 * Arguments:
74 * ==========
75 *
76 *> \param[in] LEFT
77 *> \verbatim
78 *> LEFT is LOGICAL
79 *> =.TRUE.: The eigenvectors in the columns of E are assumed
80 *> to be *left* eigenvectors.
81 *> =.FALSE.: The eigenvectors in the columns of E are assumed
82 *> to be *right* eigenvectors.
83 *> \endverbatim
84 *>
85 *> \param[in] N
86 *> \verbatim
87 *> N is INTEGER
88 *> The size of the matrices. If it is zero, DGET52 does
89 *> nothing. It must be at least zero.
90 *> \endverbatim
91 *>
92 *> \param[in] A
93 *> \verbatim
94 *> A is DOUBLE PRECISION array, dimension (LDA, N)
95 *> The matrix A.
96 *> \endverbatim
97 *>
98 *> \param[in] LDA
99 *> \verbatim
100 *> LDA is INTEGER
101 *> The leading dimension of A. It must be at least 1
102 *> and at least N.
103 *> \endverbatim
104 *>
105 *> \param[in] B
106 *> \verbatim
107 *> B is DOUBLE PRECISION array, dimension (LDB, N)
108 *> The matrix B.
109 *> \endverbatim
110 *>
111 *> \param[in] LDB
112 *> \verbatim
113 *> LDB is INTEGER
114 *> The leading dimension of B. It must be at least 1
115 *> and at least N.
116 *> \endverbatim
117 *>
118 *> \param[in] E
119 *> \verbatim
120 *> E is DOUBLE PRECISION array, dimension (LDE, N)
121 *> The matrix of eigenvectors. It must be O( 1 ). Complex
122 *> eigenvalues and eigenvectors always come in pairs, the
123 *> eigenvalue and its conjugate being stored in adjacent
124 *> elements of ALPHAR, ALPHAI, and BETA. Thus, if a(j)/b(j)
125 *> and a(j+1)/b(j+1) are a complex conjugate pair of
126 *> generalized eigenvalues, then E(,j) contains the real part
127 *> of the eigenvector and E(,j+1) contains the imaginary part.
128 *> Note that whether E(,j) is a real eigenvector or part of a
129 *> complex one is specified by whether ALPHAI(j) is zero or not.
130 *> \endverbatim
131 *>
132 *> \param[in] LDE
133 *> \verbatim
134 *> LDE is INTEGER
135 *> The leading dimension of E. It must be at least 1 and at
136 *> least N.
137 *> \endverbatim
138 *>
139 *> \param[in] ALPHAR
140 *> \verbatim
141 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
142 *> The real parts of the values a(j) as described above, which,
143 *> along with b(j), define the generalized eigenvalues.
144 *> Complex eigenvalues always come in complex conjugate pairs
145 *> a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
146 *> elements in ALPHAR, ALPHAI, and BETA. Thus, if the j-th
147 *> and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
148 *> is assumed to be equal to ALPHAR(j)/BETA(j).
149 *> \endverbatim
150 *>
151 *> \param[in] ALPHAI
152 *> \verbatim
153 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
154 *> The imaginary parts of the values a(j) as described above,
155 *> which, along with b(j), define the generalized eigenvalues.
156 *> If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
157 *> is part of a complex conjugate pair. Complex eigenvalues
158 *> always come in complex conjugate pairs a(j)/b(j) and
159 *> a(j+1)/b(j+1), which are stored in adjacent elements in
160 *> ALPHAR, ALPHAI, and BETA. Thus, if the j-th and (j+1)-st
161 *> eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
162 *> be equal to -ALPHAI(j)/BETA(j). Also, nonzero values in
163 *> ALPHAI are assumed to always come in adjacent pairs.
164 *> \endverbatim
165 *>
166 *> \param[in] BETA
167 *> \verbatim
168 *> BETA is DOUBLE PRECISION array, dimension (N)
169 *> The values b(j) as described above, which, along with a(j),
170 *> define the generalized eigenvalues.
171 *> \endverbatim
172 *>
173 *> \param[out] WORK
174 *> \verbatim
175 *> WORK is DOUBLE PRECISION array, dimension (N**2+N)
176 *> \endverbatim
177 *>
178 *> \param[out] RESULT
179 *> \verbatim
180 *> RESULT is DOUBLE PRECISION array, dimension (2)
181 *> The values computed by the test described above. If A E or
182 *> B E is likely to overflow, then RESULT(1:2) is set to
183 *> 10 / ulp.
184 *> \endverbatim
185 *
186 * Authors:
187 * ========
188 *
189 *> \author Univ. of Tennessee
190 *> \author Univ. of California Berkeley
191 *> \author Univ. of Colorado Denver
192 *> \author NAG Ltd.
193 *
194 *> \date November 2011
195 *
196 *> \ingroup double_eig
197 *
198 * =====================================================================
199  SUBROUTINE dget52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
200  $ alphai, beta, work, result )
201 *
202 * -- LAPACK test routine (version 3.4.0) --
203 * -- LAPACK is a software package provided by Univ. of Tennessee, --
204 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205 * November 2011
206 *
207 * .. Scalar Arguments ..
208  LOGICAL left
209  INTEGER lda, ldb, lde, n
210 * ..
211 * .. Array Arguments ..
212  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
213  $ b( ldb, * ), beta( * ), e( lde, * ),
214  $ result( 2 ), work( * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  DOUBLE PRECISION zero, one, ten
221  parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
222 * ..
223 * .. Local Scalars ..
224  LOGICAL ilcplx
225  CHARACTER normab, trans
226  INTEGER j, jvec
227  DOUBLE PRECISION abmax, acoef, alfmax, anorm, bcoefi, bcoefr,
228  $ betmax, bnorm, enorm, enrmer, errnrm, safmax,
229  $ safmin, salfi, salfr, sbeta, scale, temp1, ulp
230 * ..
231 * .. External Functions ..
232  DOUBLE PRECISION dlamch, dlange
233  EXTERNAL dlamch, dlange
234 * ..
235 * .. External Subroutines ..
236  EXTERNAL dgemv
237 * ..
238 * .. Intrinsic Functions ..
239  INTRINSIC abs, dble, max
240 * ..
241 * .. Executable Statements ..
242 *
243  result( 1 ) = zero
244  result( 2 ) = zero
245  IF( n.LE.0 )
246  $ return
247 *
248  safmin = dlamch( 'Safe minimum' )
249  safmax = one / safmin
250  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
251 *
252  IF( left ) THEN
253  trans = 'T'
254  normab = 'I'
255  ELSE
256  trans = 'N'
257  normab = 'O'
258  END IF
259 *
260 * Norm of A, B, and E:
261 *
262  anorm = max( dlange( normab, n, n, a, lda, work ), safmin )
263  bnorm = max( dlange( normab, n, n, b, ldb, work ), safmin )
264  enorm = max( dlange( 'O', n, n, e, lde, work ), ulp )
265  alfmax = safmax / max( one, bnorm )
266  betmax = safmax / max( one, anorm )
267 *
268 * Compute error matrix.
269 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
270 *
271  ilcplx = .false.
272  DO 10 jvec = 1, n
273  IF( ilcplx ) THEN
274 *
275 * 2nd Eigenvalue/-vector of pair -- do nothing
276 *
277  ilcplx = .false.
278  ELSE
279  salfr = alphar( jvec )
280  salfi = alphai( jvec )
281  sbeta = beta( jvec )
282  IF( salfi.EQ.zero ) THEN
283 *
284 * Real eigenvalue and -vector
285 *
286  abmax = max( abs( salfr ), abs( sbeta ) )
287  IF( abs( salfr ).GT.alfmax .OR. abs( sbeta ).GT.
288  $ betmax .OR. abmax.LT.one ) THEN
289  scale = one / max( abmax, safmin )
290  salfr = scale*salfr
291  sbeta = scale*sbeta
292  END IF
293  scale = one / max( abs( salfr )*bnorm,
294  $ abs( sbeta )*anorm, safmin )
295  acoef = scale*sbeta
296  bcoefr = scale*salfr
297  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
298  $ zero, work( n*( jvec-1 )+1 ), 1 )
299  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
300  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
301  ELSE
302 *
303 * Complex conjugate pair
304 *
305  ilcplx = .true.
306  IF( jvec.EQ.n ) THEN
307  result( 1 ) = ten / ulp
308  return
309  END IF
310  abmax = max( abs( salfr )+abs( salfi ), abs( sbeta ) )
311  IF( abs( salfr )+abs( salfi ).GT.alfmax .OR.
312  $ abs( sbeta ).GT.betmax .OR. abmax.LT.one ) THEN
313  scale = one / max( abmax, safmin )
314  salfr = scale*salfr
315  salfi = scale*salfi
316  sbeta = scale*sbeta
317  END IF
318  scale = one / max( ( abs( salfr )+abs( salfi ) )*bnorm,
319  $ abs( sbeta )*anorm, safmin )
320  acoef = scale*sbeta
321  bcoefr = scale*salfr
322  bcoefi = scale*salfi
323  IF( left ) THEN
324  bcoefi = -bcoefi
325  END IF
326 *
327  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
328  $ zero, work( n*( jvec-1 )+1 ), 1 )
329  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
330  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
331  CALL dgemv( trans, n, n, bcoefi, b, lda, e( 1, jvec+1 ),
332  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
333 *
334  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec+1 ),
335  $ 1, zero, work( n*jvec+1 ), 1 )
336  CALL dgemv( trans, n, n, -bcoefi, b, lda, e( 1, jvec ),
337  $ 1, one, work( n*jvec+1 ), 1 )
338  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec+1 ),
339  $ 1, one, work( n*jvec+1 ), 1 )
340  END IF
341  END IF
342  10 continue
343 *
344  errnrm = dlange( 'One', n, n, work, n, work( n**2+1 ) ) / enorm
345 *
346 * Compute RESULT(1)
347 *
348  result( 1 ) = errnrm / ulp
349 *
350 * Normalization of E:
351 *
352  enrmer = zero
353  ilcplx = .false.
354  DO 40 jvec = 1, n
355  IF( ilcplx ) THEN
356  ilcplx = .false.
357  ELSE
358  temp1 = zero
359  IF( alphai( jvec ).EQ.zero ) THEN
360  DO 20 j = 1, n
361  temp1 = max( temp1, abs( e( j, jvec ) ) )
362  20 continue
363  enrmer = max( enrmer, temp1-one )
364  ELSE
365  ilcplx = .true.
366  DO 30 j = 1, n
367  temp1 = max( temp1, abs( e( j, jvec ) )+
368  $ abs( e( j, jvec+1 ) ) )
369  30 continue
370  enrmer = max( enrmer, temp1-one )
371  END IF
372  END IF
373  40 continue
374 *
375 * Compute RESULT(2) : the normalization error in E.
376 *
377  result( 2 ) = enrmer / ( dble( n )*ulp )
378 *
379  return
380 *
381 * End of DGET52
382 *
383  END