LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dchkgg.f
Go to the documentation of this file.
1 *> \brief \b DCHKGG
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 DCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
13 * S2, P1, P2, U, LDU, V, Q, Z, ALPHR1, ALPHI1,
14 * BETA1, ALPHR3, ALPHI3, BETA3, EVECTL, EVECTR,
15 * WORK, LWORK, LLWORK, RESULT, INFO )
16 *
17 * .. Scalar Arguments ..
18 * LOGICAL TSTDIF
19 * INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
20 * DOUBLE PRECISION THRESH, THRSHN
21 * ..
22 * .. Array Arguments ..
23 * LOGICAL DOTYPE( * ), LLWORK( * )
24 * INTEGER ISEED( 4 ), NN( * )
25 * DOUBLE PRECISION A( LDA, * ), ALPHI1( * ), ALPHI3( * ),
26 * $ ALPHR1( * ), ALPHR3( * ), B( LDA, * ),
27 * $ BETA1( * ), BETA3( * ), EVECTL( LDU, * ),
28 * $ EVECTR( LDU, * ), H( LDA, * ), P1( LDA, * ),
29 * $ P2( LDA, * ), Q( LDU, * ), RESULT( 15 ),
30 * $ S1( LDA, * ), S2( LDA, * ), T( LDA, * ),
31 * $ U( LDU, * ), V( LDU, * ), WORK( * ),
32 * $ Z( LDU, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DCHKGG checks the nonsymmetric generalized eigenvalue problem
42 *> routines.
43 *> T T T
44 *> DGGHRD factors A and B as U H V and U T V , where means
45 *> transpose, H is hessenberg, T is triangular and U and V are
46 *> orthogonal.
47 *> T T
48 *> DHGEQZ factors H and T as Q S Z and Q P Z , where P is upper
49 *> triangular, S is in generalized Schur form (block upper triangular,
50 *> with 1x1 and 2x2 blocks on the diagonal, the 2x2 blocks
51 *> corresponding to complex conjugate pairs of generalized
52 *> eigenvalues), and Q and Z are orthogonal. It also computes the
53 *> generalized eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)),
54 *> where alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus,
55 *> w(j) = alpha(j)/beta(j) is a root of the generalized eigenvalue
56 *> problem
57 *>
58 *> det( A - w(j) B ) = 0
59 *>
60 *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
61 *> problem
62 *>
63 *> det( m(j) A - B ) = 0
64 *>
65 *> DTGEVC computes the matrix L of left eigenvectors and the matrix R
66 *> of right eigenvectors for the matrix pair ( S, P ). In the
67 *> description below, l and r are left and right eigenvectors
68 *> corresponding to the generalized eigenvalues (alpha,beta).
69 *>
70 *> When DCHKGG is called, a number of matrix "sizes" ("n's") and a
71 *> number of matrix "types" are specified. For each size ("n")
72 *> and each type of matrix, one matrix will be generated and used
73 *> to test the nonsymmetric eigenroutines. For each matrix, 15
74 *> tests will be performed. The first twelve "test ratios" should be
75 *> small -- O(1). They will be compared with the threshhold THRESH:
76 *>
77 *> T
78 *> (1) | A - U H V | / ( |A| n ulp )
79 *>
80 *> T
81 *> (2) | B - U T V | / ( |B| n ulp )
82 *>
83 *> T
84 *> (3) | I - UU | / ( n ulp )
85 *>
86 *> T
87 *> (4) | I - VV | / ( n ulp )
88 *>
89 *> T
90 *> (5) | H - Q S Z | / ( |H| n ulp )
91 *>
92 *> T
93 *> (6) | T - Q P Z | / ( |T| n ulp )
94 *>
95 *> T
96 *> (7) | I - QQ | / ( n ulp )
97 *>
98 *> T
99 *> (8) | I - ZZ | / ( n ulp )
100 *>
101 *> (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
102 *>
103 *> | l**H * (beta S - alpha P) | / ( ulp max( |beta S|, |alpha P| ) )
104 *>
105 *> (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of
106 *> T
107 *> | l'**H * (beta H - alpha T) | / ( ulp max( |beta H|, |alpha T| ) )
108 *>
109 *> where the eigenvectors l' are the result of passing Q to
110 *> DTGEVC and back transforming (HOWMNY='B').
111 *>
112 *> (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
113 *>
114 *> | (beta S - alpha T) r | / ( ulp max( |beta S|, |alpha T| ) )
115 *>
116 *> (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of
117 *>
118 *> | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
119 *>
120 *> where the eigenvectors r' are the result of passing Z to
121 *> DTGEVC and back transforming (HOWMNY='B').
122 *>
123 *> The last three test ratios will usually be small, but there is no
124 *> mathematical requirement that they be so. They are therefore
125 *> compared with THRESH only if TSTDIF is .TRUE.
126 *>
127 *> (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
128 *>
129 *> (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
130 *>
131 *> (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
132 *> |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
133 *>
134 *> In addition, the normalization of L and R are checked, and compared
135 *> with the threshhold THRSHN.
136 *>
137 *> Test Matrices
138 *> ---- --------
139 *>
140 *> The sizes of the test matrices are specified by an array
141 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
142 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
143 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
144 *> Currently, the list of possible types is:
145 *>
146 *> (1) ( 0, 0 ) (a pair of zero matrices)
147 *>
148 *> (2) ( I, 0 ) (an identity and a zero matrix)
149 *>
150 *> (3) ( 0, I ) (an identity and a zero matrix)
151 *>
152 *> (4) ( I, I ) (a pair of identity matrices)
153 *>
154 *> t t
155 *> (5) ( J , J ) (a pair of transposed Jordan blocks)
156 *>
157 *> t ( I 0 )
158 *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
159 *> ( 0 I ) ( 0 J )
160 *> and I is a k x k identity and J a (k+1)x(k+1)
161 *> Jordan block; k=(N-1)/2
162 *>
163 *> (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
164 *> matrix with those diagonal entries.)
165 *> (8) ( I, D )
166 *>
167 *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
168 *>
169 *> (10) ( small*D, big*I )
170 *>
171 *> (11) ( big*I, small*D )
172 *>
173 *> (12) ( small*I, big*D )
174 *>
175 *> (13) ( big*D, big*I )
176 *>
177 *> (14) ( small*D, small*I )
178 *>
179 *> (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
180 *> D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
181 *> t t
182 *> (16) U ( J , J ) V where U and V are random orthogonal matrices.
183 *>
184 *> (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices
185 *> with random O(1) entries above the diagonal
186 *> and diagonal entries diag(T1) =
187 *> ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
188 *> ( 0, N-3, N-4,..., 1, 0, 0 )
189 *>
190 *> (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
191 *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
192 *> s = machine precision.
193 *>
194 *> (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
195 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
196 *>
197 *> N-5
198 *> (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
199 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
200 *>
201 *> (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
202 *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
203 *> where r1,..., r(N-4) are random.
204 *>
205 *> (22) U ( big*T1, small*T2 ) V diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
206 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
207 *>
208 *> (23) U ( small*T1, big*T2 ) V diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
209 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
210 *>
211 *> (24) U ( small*T1, small*T2 ) V diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
212 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
213 *>
214 *> (25) U ( big*T1, big*T2 ) V diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
215 *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
216 *>
217 *> (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular
218 *> matrices.
219 *> \endverbatim
220 *
221 * Arguments:
222 * ==========
223 *
224 *> \param[in] NSIZES
225 *> \verbatim
226 *> NSIZES is INTEGER
227 *> The number of sizes of matrices to use. If it is zero,
228 *> DCHKGG does nothing. It must be at least zero.
229 *> \endverbatim
230 *>
231 *> \param[in] NN
232 *> \verbatim
233 *> NN is INTEGER array, dimension (NSIZES)
234 *> An array containing the sizes to be used for the matrices.
235 *> Zero values will be skipped. The values must be at least
236 *> zero.
237 *> \endverbatim
238 *>
239 *> \param[in] NTYPES
240 *> \verbatim
241 *> NTYPES is INTEGER
242 *> The number of elements in DOTYPE. If it is zero, DCHKGG
243 *> does nothing. It must be at least zero. If it is MAXTYP+1
244 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
245 *> defined, which is to use whatever matrix is in A. This
246 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
247 *> DOTYPE(MAXTYP+1) is .TRUE. .
248 *> \endverbatim
249 *>
250 *> \param[in] DOTYPE
251 *> \verbatim
252 *> DOTYPE is LOGICAL array, dimension (NTYPES)
253 *> If DOTYPE(j) is .TRUE., then for each size in NN a
254 *> matrix of that size and of type j will be generated.
255 *> If NTYPES is smaller than the maximum number of types
256 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
257 *> MAXTYP will not be generated. If NTYPES is larger
258 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
259 *> will be ignored.
260 *> \endverbatim
261 *>
262 *> \param[in,out] ISEED
263 *> \verbatim
264 *> ISEED is INTEGER array, dimension (4)
265 *> On entry ISEED specifies the seed of the random number
266 *> generator. The array elements should be between 0 and 4095;
267 *> if not they will be reduced mod 4096. Also, ISEED(4) must
268 *> be odd. The random number generator uses a linear
269 *> congruential sequence limited to small integers, and so
270 *> should produce machine independent random numbers. The
271 *> values of ISEED are changed on exit, and can be used in the
272 *> next call to DCHKGG to continue the same random number
273 *> sequence.
274 *> \endverbatim
275 *>
276 *> \param[in] THRESH
277 *> \verbatim
278 *> THRESH is DOUBLE PRECISION
279 *> A test will count as "failed" if the "error", computed as
280 *> described above, exceeds THRESH. Note that the error is
281 *> scaled to be O(1), so THRESH should be a reasonably small
282 *> multiple of 1, e.g., 10 or 100. In particular, it should
283 *> not depend on the precision (single vs. double) or the size
284 *> of the matrix. It must be at least zero.
285 *> \endverbatim
286 *>
287 *> \param[in] TSTDIF
288 *> \verbatim
289 *> TSTDIF is LOGICAL
290 *> Specifies whether test ratios 13-15 will be computed and
291 *> compared with THRESH.
292 *> = .FALSE.: Only test ratios 1-12 will be computed and tested.
293 *> Ratios 13-15 will be set to zero.
294 *> = .TRUE.: All the test ratios 1-15 will be computed and
295 *> tested.
296 *> \endverbatim
297 *>
298 *> \param[in] THRSHN
299 *> \verbatim
300 *> THRSHN is DOUBLE PRECISION
301 *> Threshhold for reporting eigenvector normalization error.
302 *> If the normalization of any eigenvector differs from 1 by
303 *> more than THRSHN*ulp, then a special error message will be
304 *> printed. (This is handled separately from the other tests,
305 *> since only a compiler or programming error should cause an
306 *> error message, at least if THRSHN is at least 5--10.)
307 *> \endverbatim
308 *>
309 *> \param[in] NOUNIT
310 *> \verbatim
311 *> NOUNIT is INTEGER
312 *> The FORTRAN unit number for printing out error messages
313 *> (e.g., if a routine returns IINFO not equal to 0.)
314 *> \endverbatim
315 *>
316 *> \param[in,out] A
317 *> \verbatim
318 *> A is DOUBLE PRECISION array, dimension
319 *> (LDA, max(NN))
320 *> Used to hold the original A matrix. Used as input only
321 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
322 *> DOTYPE(MAXTYP+1)=.TRUE.
323 *> \endverbatim
324 *>
325 *> \param[in] LDA
326 *> \verbatim
327 *> LDA is INTEGER
328 *> The leading dimension of A, B, H, T, S1, P1, S2, and P2.
329 *> It must be at least 1 and at least max( NN ).
330 *> \endverbatim
331 *>
332 *> \param[in,out] B
333 *> \verbatim
334 *> B is DOUBLE PRECISION array, dimension
335 *> (LDA, max(NN))
336 *> Used to hold the original B matrix. Used as input only
337 *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
338 *> DOTYPE(MAXTYP+1)=.TRUE.
339 *> \endverbatim
340 *>
341 *> \param[out] H
342 *> \verbatim
343 *> H is DOUBLE PRECISION array, dimension (LDA, max(NN))
344 *> The upper Hessenberg matrix computed from A by DGGHRD.
345 *> \endverbatim
346 *>
347 *> \param[out] T
348 *> \verbatim
349 *> T is DOUBLE PRECISION array, dimension (LDA, max(NN))
350 *> The upper triangular matrix computed from B by DGGHRD.
351 *> \endverbatim
352 *>
353 *> \param[out] S1
354 *> \verbatim
355 *> S1 is DOUBLE PRECISION array, dimension (LDA, max(NN))
356 *> The Schur (block upper triangular) matrix computed from H by
357 *> DHGEQZ when Q and Z are also computed.
358 *> \endverbatim
359 *>
360 *> \param[out] S2
361 *> \verbatim
362 *> S2 is DOUBLE PRECISION array, dimension (LDA, max(NN))
363 *> The Schur (block upper triangular) matrix computed from H by
364 *> DHGEQZ when Q and Z are not computed.
365 *> \endverbatim
366 *>
367 *> \param[out] P1
368 *> \verbatim
369 *> P1 is DOUBLE PRECISION array, dimension (LDA, max(NN))
370 *> The upper triangular matrix computed from T by DHGEQZ
371 *> when Q and Z are also computed.
372 *> \endverbatim
373 *>
374 *> \param[out] P2
375 *> \verbatim
376 *> P2 is DOUBLE PRECISION array, dimension (LDA, max(NN))
377 *> The upper triangular matrix computed from T by DHGEQZ
378 *> when Q and Z are not computed.
379 *> \endverbatim
380 *>
381 *> \param[out] U
382 *> \verbatim
383 *> U is DOUBLE PRECISION array, dimension (LDU, max(NN))
384 *> The (left) orthogonal matrix computed by DGGHRD.
385 *> \endverbatim
386 *>
387 *> \param[in] LDU
388 *> \verbatim
389 *> LDU is INTEGER
390 *> The leading dimension of U, V, Q, Z, EVECTL, and EVEZTR. It
391 *> must be at least 1 and at least max( NN ).
392 *> \endverbatim
393 *>
394 *> \param[out] V
395 *> \verbatim
396 *> V is DOUBLE PRECISION array, dimension (LDU, max(NN))
397 *> The (right) orthogonal matrix computed by DGGHRD.
398 *> \endverbatim
399 *>
400 *> \param[out] Q
401 *> \verbatim
402 *> Q is DOUBLE PRECISION array, dimension (LDU, max(NN))
403 *> The (left) orthogonal matrix computed by DHGEQZ.
404 *> \endverbatim
405 *>
406 *> \param[out] Z
407 *> \verbatim
408 *> Z is DOUBLE PRECISION array, dimension (LDU, max(NN))
409 *> The (left) orthogonal matrix computed by DHGEQZ.
410 *> \endverbatim
411 *>
412 *> \param[out] ALPHR1
413 *> \verbatim
414 *> ALPHR1 is DOUBLE PRECISION array, dimension (max(NN))
415 *> \endverbatim
416 *>
417 *> \param[out] ALPHI1
418 *> \verbatim
419 *> ALPHI1 is DOUBLE PRECISION array, dimension (max(NN))
420 *> \endverbatim
421 *>
422 *> \param[out] BETA1
423 *> \verbatim
424 *> BETA1 is DOUBLE PRECISION array, dimension (max(NN))
425 *>
426 *> The generalized eigenvalues of (A,B) computed by DHGEQZ
427 *> when Q, Z, and the full Schur matrices are computed.
428 *> On exit, ( ALPHR1(k)+ALPHI1(k)*i ) / BETA1(k) is the k-th
429 *> generalized eigenvalue of the matrices in A and B.
430 *> \endverbatim
431 *>
432 *> \param[out] ALPHR3
433 *> \verbatim
434 *> ALPHR3 is DOUBLE PRECISION array, dimension (max(NN))
435 *> \endverbatim
436 *>
437 *> \param[out] ALPHI3
438 *> \verbatim
439 *> ALPHI3 is DOUBLE PRECISION array, dimension (max(NN))
440 *> \endverbatim
441 *>
442 *> \param[out] BETA3
443 *> \verbatim
444 *> BETA3 is DOUBLE PRECISION array, dimension (max(NN))
445 *> \endverbatim
446 *>
447 *> \param[out] EVECTL
448 *> \verbatim
449 *> EVECTL is DOUBLE PRECISION array, dimension (LDU, max(NN))
450 *> The (block lower triangular) left eigenvector matrix for
451 *> the matrices in S1 and P1. (See DTGEVC for the format.)
452 *> \endverbatim
453 *>
454 *> \param[out] EVECTR
455 *> \verbatim
456 *> EVECTR is DOUBLE PRECISION array, dimension (LDU, max(NN))
457 *> The (block upper triangular) right eigenvector matrix for
458 *> the matrices in S1 and P1. (See DTGEVC for the format.)
459 *> \endverbatim
460 *>
461 *> \param[out] WORK
462 *> \verbatim
463 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
464 *> \endverbatim
465 *>
466 *> \param[in] LWORK
467 *> \verbatim
468 *> LWORK is INTEGER
469 *> The number of entries in WORK. This must be at least
470 *> max( 2 * N**2, 6*N, 1 ), for all N=NN(j).
471 *> \endverbatim
472 *>
473 *> \param[out] LLWORK
474 *> \verbatim
475 *> LLWORK is LOGICAL array, dimension (max(NN))
476 *> \endverbatim
477 *>
478 *> \param[out] RESULT
479 *> \verbatim
480 *> RESULT is DOUBLE PRECISION array, dimension (15)
481 *> The values computed by the tests described above.
482 *> The values are currently limited to 1/ulp, to avoid
483 *> overflow.
484 *> \endverbatim
485 *>
486 *> \param[out] INFO
487 *> \verbatim
488 *> INFO is INTEGER
489 *> = 0: successful exit
490 *> < 0: if INFO = -i, the i-th argument had an illegal value
491 *> > 0: A routine returned an error code. INFO is the
492 *> absolute value of the INFO value returned.
493 *> \endverbatim
494 *
495 * Authors:
496 * ========
497 *
498 *> \author Univ. of Tennessee
499 *> \author Univ. of California Berkeley
500 *> \author Univ. of Colorado Denver
501 *> \author NAG Ltd.
502 *
503 *> \date November 2011
504 *
505 *> \ingroup double_eig
506 *
507 * =====================================================================
508  SUBROUTINE dchkgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
509  $ tstdif, thrshn, nounit, a, lda, b, h, t, s1,
510  $ s2, p1, p2, u, ldu, v, q, z, alphr1, alphi1,
511  $ beta1, alphr3, alphi3, beta3, evectl, evectr,
512  $ work, lwork, llwork, result, info )
513 *
514 * -- LAPACK test routine (version 3.4.0) --
515 * -- LAPACK is a software package provided by Univ. of Tennessee, --
516 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
517 * November 2011
518 *
519 * .. Scalar Arguments ..
520  LOGICAL tstdif
521  INTEGER info, lda, ldu, lwork, nounit, nsizes, ntypes
522  DOUBLE PRECISION thresh, thrshn
523 * ..
524 * .. Array Arguments ..
525  LOGICAL dotype( * ), llwork( * )
526  INTEGER iseed( 4 ), nn( * )
527  DOUBLE PRECISION a( lda, * ), alphi1( * ), alphi3( * ),
528  $ alphr1( * ), alphr3( * ), b( lda, * ),
529  $ beta1( * ), beta3( * ), evectl( ldu, * ),
530  $ evectr( ldu, * ), h( lda, * ), p1( lda, * ),
531  $ p2( lda, * ), q( ldu, * ), result( 15 ),
532  $ s1( lda, * ), s2( lda, * ), t( lda, * ),
533  $ u( ldu, * ), v( ldu, * ), work( * ),
534  $ z( ldu, * )
535 * ..
536 *
537 * =====================================================================
538 *
539 * .. Parameters ..
540  DOUBLE PRECISION zero, one
541  parameter( zero = 0.0d0, one = 1.0d0 )
542  INTEGER maxtyp
543  parameter( maxtyp = 26 )
544 * ..
545 * .. Local Scalars ..
546  LOGICAL badnn
547  INTEGER i1, iadd, iinfo, in, j, jc, jr, jsize, jtype,
548  $ lwkopt, mtypes, n, n1, nerrs, nmats, nmax,
549  $ ntest, ntestt
550  DOUBLE PRECISION anorm, bnorm, safmax, safmin, temp1, temp2,
551  $ ulp, ulpinv
552 * ..
553 * .. Local Arrays ..
554  INTEGER iasign( maxtyp ), ibsign( maxtyp ),
555  $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
556  $ katype( maxtyp ), kazero( maxtyp ),
557  $ kbmagn( maxtyp ), kbtype( maxtyp ),
558  $ kbzero( maxtyp ), kclass( maxtyp ),
559  $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
560  DOUBLE PRECISION dumma( 4 ), rmagn( 0: 3 )
561 * ..
562 * .. External Functions ..
563  DOUBLE PRECISION dlamch, dlange, dlarnd
564  EXTERNAL dlamch, dlange, dlarnd
565 * ..
566 * .. External Subroutines ..
567  EXTERNAL dgeqr2, dget51, dget52, dgghrd, dhgeqz, dlabad,
569  $ dtgevc, xerbla
570 * ..
571 * .. Intrinsic Functions ..
572  INTRINSIC abs, dble, max, min, sign
573 * ..
574 * .. Data statements ..
575  DATA kclass / 15*1, 10*2, 1*3 /
576  DATA kz1 / 0, 1, 2, 1, 3, 3 /
577  DATA kz2 / 0, 0, 1, 2, 1, 1 /
578  DATA kadd / 0, 0, 0, 0, 3, 2 /
579  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
580  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
581  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
582  $ 1, 1, -4, 2, -4, 8*8, 0 /
583  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
584  $ 4*5, 4*3, 1 /
585  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
586  $ 4*6, 4*4, 1 /
587  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
588  $ 2, 1 /
589  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
590  $ 2, 1 /
591  DATA ktrian / 16*0, 10*1 /
592  DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
593  $ 5*2, 0 /
594  DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
595 * ..
596 * .. Executable Statements ..
597 *
598 * Check for errors
599 *
600  info = 0
601 *
602  badnn = .false.
603  nmax = 1
604  DO 10 j = 1, nsizes
605  nmax = max( nmax, nn( j ) )
606  IF( nn( j ).LT.0 )
607  $ badnn = .true.
608  10 continue
609 *
610 * Maximum blocksize and shift -- we assume that blocksize and number
611 * of shifts are monotone increasing functions of N.
612 *
613  lwkopt = max( 6*nmax, 2*nmax*nmax, 1 )
614 *
615 * Check for errors
616 *
617  IF( nsizes.LT.0 ) THEN
618  info = -1
619  ELSE IF( badnn ) THEN
620  info = -2
621  ELSE IF( ntypes.LT.0 ) THEN
622  info = -3
623  ELSE IF( thresh.LT.zero ) THEN
624  info = -6
625  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
626  info = -10
627  ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax ) THEN
628  info = -19
629  ELSE IF( lwkopt.GT.lwork ) THEN
630  info = -30
631  END IF
632 *
633  IF( info.NE.0 ) THEN
634  CALL xerbla( 'DCHKGG', -info )
635  return
636  END IF
637 *
638 * Quick return if possible
639 *
640  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
641  $ return
642 *
643  safmin = dlamch( 'Safe minimum' )
644  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
645  safmin = safmin / ulp
646  safmax = one / safmin
647  CALL dlabad( safmin, safmax )
648  ulpinv = one / ulp
649 *
650 * The values RMAGN(2:3) depend on N, see below.
651 *
652  rmagn( 0 ) = zero
653  rmagn( 1 ) = one
654 *
655 * Loop over sizes, types
656 *
657  ntestt = 0
658  nerrs = 0
659  nmats = 0
660 *
661  DO 240 jsize = 1, nsizes
662  n = nn( jsize )
663  n1 = max( 1, n )
664  rmagn( 2 ) = safmax*ulp / dble( n1 )
665  rmagn( 3 ) = safmin*ulpinv*n1
666 *
667  IF( nsizes.NE.1 ) THEN
668  mtypes = min( maxtyp, ntypes )
669  ELSE
670  mtypes = min( maxtyp+1, ntypes )
671  END IF
672 *
673  DO 230 jtype = 1, mtypes
674  IF( .NOT.dotype( jtype ) )
675  $ go to 230
676  nmats = nmats + 1
677  ntest = 0
678 *
679 * Save ISEED in case of an error.
680 *
681  DO 20 j = 1, 4
682  ioldsd( j ) = iseed( j )
683  20 continue
684 *
685 * Initialize RESULT
686 *
687  DO 30 j = 1, 15
688  result( j ) = zero
689  30 continue
690 *
691 * Compute A and B
692 *
693 * Description of control parameters:
694 *
695 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
696 * =3 means random.
697 * KATYPE: the "type" to be passed to DLATM4 for computing A.
698 * KAZERO: the pattern of zeros on the diagonal for A:
699 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
700 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
701 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
702 * non-zero entries.)
703 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
704 * =2: large, =3: small.
705 * IASIGN: 1 if the diagonal elements of A are to be
706 * multiplied by a random magnitude 1 number, =2 if
707 * randomly chosen diagonal blocks are to be rotated
708 * to form 2x2 blocks.
709 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
710 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
711 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
712 * RMAGN: used to implement KAMAGN and KBMAGN.
713 *
714  IF( mtypes.GT.maxtyp )
715  $ go to 110
716  iinfo = 0
717  IF( kclass( jtype ).LT.3 ) THEN
718 *
719 * Generate A (w/o rotation)
720 *
721  IF( abs( katype( jtype ) ).EQ.3 ) THEN
722  in = 2*( ( n-1 ) / 2 ) + 1
723  IF( in.NE.n )
724  $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
725  ELSE
726  in = n
727  END IF
728  CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
729  $ kz2( kazero( jtype ) ), iasign( jtype ),
730  $ rmagn( kamagn( jtype ) ), ulp,
731  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
732  $ iseed, a, lda )
733  iadd = kadd( kazero( jtype ) )
734  IF( iadd.GT.0 .AND. iadd.LE.n )
735  $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
736 *
737 * Generate B (w/o rotation)
738 *
739  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
740  in = 2*( ( n-1 ) / 2 ) + 1
741  IF( in.NE.n )
742  $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
743  ELSE
744  in = n
745  END IF
746  CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
747  $ kz2( kbzero( jtype ) ), ibsign( jtype ),
748  $ rmagn( kbmagn( jtype ) ), one,
749  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
750  $ iseed, b, lda )
751  iadd = kadd( kbzero( jtype ) )
752  IF( iadd.NE.0 .AND. iadd.LE.n )
753  $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
754 *
755  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
756 *
757 * Include rotations
758 *
759 * Generate U, V as Householder transformations times
760 * a diagonal matrix.
761 *
762  DO 50 jc = 1, n - 1
763  DO 40 jr = jc, n
764  u( jr, jc ) = dlarnd( 3, iseed )
765  v( jr, jc ) = dlarnd( 3, iseed )
766  40 continue
767  CALL dlarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
768  $ work( jc ) )
769  work( 2*n+jc ) = sign( one, u( jc, jc ) )
770  u( jc, jc ) = one
771  CALL dlarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
772  $ work( n+jc ) )
773  work( 3*n+jc ) = sign( one, v( jc, jc ) )
774  v( jc, jc ) = one
775  50 continue
776  u( n, n ) = one
777  work( n ) = zero
778  work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
779  v( n, n ) = one
780  work( 2*n ) = zero
781  work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
782 *
783 * Apply the diagonal matrices
784 *
785  DO 70 jc = 1, n
786  DO 60 jr = 1, n
787  a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
788  $ a( jr, jc )
789  b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
790  $ b( jr, jc )
791  60 continue
792  70 continue
793  CALL dorm2r( 'L', 'N', n, n, n-1, u, ldu, work, a,
794  $ lda, work( 2*n+1 ), iinfo )
795  IF( iinfo.NE.0 )
796  $ go to 100
797  CALL dorm2r( 'R', 'T', n, n, n-1, v, ldu, work( n+1 ),
798  $ a, lda, work( 2*n+1 ), iinfo )
799  IF( iinfo.NE.0 )
800  $ go to 100
801  CALL dorm2r( 'L', 'N', n, n, n-1, u, ldu, work, b,
802  $ lda, work( 2*n+1 ), iinfo )
803  IF( iinfo.NE.0 )
804  $ go to 100
805  CALL dorm2r( 'R', 'T', n, n, n-1, v, ldu, work( n+1 ),
806  $ b, lda, work( 2*n+1 ), iinfo )
807  IF( iinfo.NE.0 )
808  $ go to 100
809  END IF
810  ELSE
811 *
812 * Random matrices
813 *
814  DO 90 jc = 1, n
815  DO 80 jr = 1, n
816  a( jr, jc ) = rmagn( kamagn( jtype ) )*
817  $ dlarnd( 2, iseed )
818  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
819  $ dlarnd( 2, iseed )
820  80 continue
821  90 continue
822  END IF
823 *
824  anorm = dlange( '1', n, n, a, lda, work )
825  bnorm = dlange( '1', n, n, b, lda, work )
826 *
827  100 continue
828 *
829  IF( iinfo.NE.0 ) THEN
830  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
831  $ ioldsd
832  info = abs( iinfo )
833  return
834  END IF
835 *
836  110 continue
837 *
838 * Call DGEQR2, DORM2R, and DGGHRD to compute H, T, U, and V
839 *
840  CALL dlacpy( ' ', n, n, a, lda, h, lda )
841  CALL dlacpy( ' ', n, n, b, lda, t, lda )
842  ntest = 1
843  result( 1 ) = ulpinv
844 *
845  CALL dgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
846  IF( iinfo.NE.0 ) THEN
847  WRITE( nounit, fmt = 9999 )'DGEQR2', iinfo, n, jtype,
848  $ ioldsd
849  info = abs( iinfo )
850  go to 210
851  END IF
852 *
853  CALL dorm2r( 'L', 'T', n, n, n, t, lda, work, h, lda,
854  $ work( n+1 ), iinfo )
855  IF( iinfo.NE.0 ) THEN
856  WRITE( nounit, fmt = 9999 )'DORM2R', iinfo, n, jtype,
857  $ ioldsd
858  info = abs( iinfo )
859  go to 210
860  END IF
861 *
862  CALL dlaset( 'Full', n, n, zero, one, u, ldu )
863  CALL dorm2r( 'R', 'N', n, n, n, t, lda, work, u, ldu,
864  $ work( n+1 ), iinfo )
865  IF( iinfo.NE.0 ) THEN
866  WRITE( nounit, fmt = 9999 )'DORM2R', iinfo, n, jtype,
867  $ ioldsd
868  info = abs( iinfo )
869  go to 210
870  END IF
871 *
872  CALL dgghrd( 'V', 'I', n, 1, n, h, lda, t, lda, u, ldu, v,
873  $ ldu, iinfo )
874  IF( iinfo.NE.0 ) THEN
875  WRITE( nounit, fmt = 9999 )'DGGHRD', iinfo, n, jtype,
876  $ ioldsd
877  info = abs( iinfo )
878  go to 210
879  END IF
880  ntest = 4
881 *
882 * Do tests 1--4
883 *
884  CALL dget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
885  $ result( 1 ) )
886  CALL dget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
887  $ result( 2 ) )
888  CALL dget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
889  $ result( 3 ) )
890  CALL dget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
891  $ result( 4 ) )
892 *
893 * Call DHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
894 *
895 * Compute T1 and UZ
896 *
897 * Eigenvalues only
898 *
899  CALL dlacpy( ' ', n, n, h, lda, s2, lda )
900  CALL dlacpy( ' ', n, n, t, lda, p2, lda )
901  ntest = 5
902  result( 5 ) = ulpinv
903 *
904  CALL dhgeqz( 'E', 'N', 'N', n, 1, n, s2, lda, p2, lda,
905  $ alphr3, alphi3, beta3, q, ldu, z, ldu, work,
906  $ lwork, iinfo )
907  IF( iinfo.NE.0 ) THEN
908  WRITE( nounit, fmt = 9999 )'DHGEQZ(E)', iinfo, n, jtype,
909  $ ioldsd
910  info = abs( iinfo )
911  go to 210
912  END IF
913 *
914 * Eigenvalues and Full Schur Form
915 *
916  CALL dlacpy( ' ', n, n, h, lda, s2, lda )
917  CALL dlacpy( ' ', n, n, t, lda, p2, lda )
918 *
919  CALL dhgeqz( 'S', 'N', 'N', n, 1, n, s2, lda, p2, lda,
920  $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
921  $ lwork, iinfo )
922  IF( iinfo.NE.0 ) THEN
923  WRITE( nounit, fmt = 9999 )'DHGEQZ(S)', iinfo, n, jtype,
924  $ ioldsd
925  info = abs( iinfo )
926  go to 210
927  END IF
928 *
929 * Eigenvalues, Schur Form, and Schur Vectors
930 *
931  CALL dlacpy( ' ', n, n, h, lda, s1, lda )
932  CALL dlacpy( ' ', n, n, t, lda, p1, lda )
933 *
934  CALL dhgeqz( 'S', 'I', 'I', n, 1, n, s1, lda, p1, lda,
935  $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
936  $ lwork, iinfo )
937  IF( iinfo.NE.0 ) THEN
938  WRITE( nounit, fmt = 9999 )'DHGEQZ(V)', iinfo, n, jtype,
939  $ ioldsd
940  info = abs( iinfo )
941  go to 210
942  END IF
943 *
944  ntest = 8
945 *
946 * Do Tests 5--8
947 *
948  CALL dget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
949  $ result( 5 ) )
950  CALL dget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
951  $ result( 6 ) )
952  CALL dget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
953  $ result( 7 ) )
954  CALL dget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
955  $ result( 8 ) )
956 *
957 * Compute the Left and Right Eigenvectors of (S1,P1)
958 *
959 * 9: Compute the left eigenvector Matrix without
960 * back transforming:
961 *
962  ntest = 9
963  result( 9 ) = ulpinv
964 *
965 * To test "SELECT" option, compute half of the eigenvectors
966 * in one call, and half in another
967 *
968  i1 = n / 2
969  DO 120 j = 1, i1
970  llwork( j ) = .true.
971  120 continue
972  DO 130 j = i1 + 1, n
973  llwork( j ) = .false.
974  130 continue
975 *
976  CALL dtgevc( 'L', 'S', llwork, n, s1, lda, p1, lda, evectl,
977  $ ldu, dumma, ldu, n, in, work, iinfo )
978  IF( iinfo.NE.0 ) THEN
979  WRITE( nounit, fmt = 9999 )'DTGEVC(L,S1)', iinfo, n,
980  $ jtype, ioldsd
981  info = abs( iinfo )
982  go to 210
983  END IF
984 *
985  i1 = in
986  DO 140 j = 1, i1
987  llwork( j ) = .false.
988  140 continue
989  DO 150 j = i1 + 1, n
990  llwork( j ) = .true.
991  150 continue
992 *
993  CALL dtgevc( 'L', 'S', llwork, n, s1, lda, p1, lda,
994  $ evectl( 1, i1+1 ), ldu, dumma, ldu, n, in,
995  $ work, iinfo )
996  IF( iinfo.NE.0 ) THEN
997  WRITE( nounit, fmt = 9999 )'DTGEVC(L,S2)', iinfo, n,
998  $ jtype, ioldsd
999  info = abs( iinfo )
1000  go to 210
1001  END IF
1002 *
1003  CALL dget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1004  $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1005  result( 9 ) = dumma( 1 )
1006  IF( dumma( 2 ).GT.thrshn ) THEN
1007  WRITE( nounit, fmt = 9998 )'Left', 'DTGEVC(HOWMNY=S)',
1008  $ dumma( 2 ), n, jtype, ioldsd
1009  END IF
1010 *
1011 * 10: Compute the left eigenvector Matrix with
1012 * back transforming:
1013 *
1014  ntest = 10
1015  result( 10 ) = ulpinv
1016  CALL dlacpy( 'F', n, n, q, ldu, evectl, ldu )
1017  CALL dtgevc( 'L', 'B', llwork, n, s1, lda, p1, lda, evectl,
1018  $ ldu, dumma, ldu, n, in, work, iinfo )
1019  IF( iinfo.NE.0 ) THEN
1020  WRITE( nounit, fmt = 9999 )'DTGEVC(L,B)', iinfo, n,
1021  $ jtype, ioldsd
1022  info = abs( iinfo )
1023  go to 210
1024  END IF
1025 *
1026  CALL dget52( .true., n, h, lda, t, lda, evectl, ldu, alphr1,
1027  $ alphi1, beta1, work, dumma( 1 ) )
1028  result( 10 ) = dumma( 1 )
1029  IF( dumma( 2 ).GT.thrshn ) THEN
1030  WRITE( nounit, fmt = 9998 )'Left', 'DTGEVC(HOWMNY=B)',
1031  $ dumma( 2 ), n, jtype, ioldsd
1032  END IF
1033 *
1034 * 11: Compute the right eigenvector Matrix without
1035 * back transforming:
1036 *
1037  ntest = 11
1038  result( 11 ) = ulpinv
1039 *
1040 * To test "SELECT" option, compute half of the eigenvectors
1041 * in one call, and half in another
1042 *
1043  i1 = n / 2
1044  DO 160 j = 1, i1
1045  llwork( j ) = .true.
1046  160 continue
1047  DO 170 j = i1 + 1, n
1048  llwork( j ) = .false.
1049  170 continue
1050 *
1051  CALL dtgevc( 'R', 'S', llwork, n, s1, lda, p1, lda, dumma,
1052  $ ldu, evectr, ldu, n, in, work, iinfo )
1053  IF( iinfo.NE.0 ) THEN
1054  WRITE( nounit, fmt = 9999 )'DTGEVC(R,S1)', iinfo, n,
1055  $ jtype, ioldsd
1056  info = abs( iinfo )
1057  go to 210
1058  END IF
1059 *
1060  i1 = in
1061  DO 180 j = 1, i1
1062  llwork( j ) = .false.
1063  180 continue
1064  DO 190 j = i1 + 1, n
1065  llwork( j ) = .true.
1066  190 continue
1067 *
1068  CALL dtgevc( 'R', 'S', llwork, n, s1, lda, p1, lda, dumma,
1069  $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1070  $ iinfo )
1071  IF( iinfo.NE.0 ) THEN
1072  WRITE( nounit, fmt = 9999 )'DTGEVC(R,S2)', iinfo, n,
1073  $ jtype, ioldsd
1074  info = abs( iinfo )
1075  go to 210
1076  END IF
1077 *
1078  CALL dget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1079  $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1080  result( 11 ) = dumma( 1 )
1081  IF( dumma( 2 ).GT.thresh ) THEN
1082  WRITE( nounit, fmt = 9998 )'Right', 'DTGEVC(HOWMNY=S)',
1083  $ dumma( 2 ), n, jtype, ioldsd
1084  END IF
1085 *
1086 * 12: Compute the right eigenvector Matrix with
1087 * back transforming:
1088 *
1089  ntest = 12
1090  result( 12 ) = ulpinv
1091  CALL dlacpy( 'F', n, n, z, ldu, evectr, ldu )
1092  CALL dtgevc( 'R', 'B', llwork, n, s1, lda, p1, lda, dumma,
1093  $ ldu, evectr, ldu, n, in, work, iinfo )
1094  IF( iinfo.NE.0 ) THEN
1095  WRITE( nounit, fmt = 9999 )'DTGEVC(R,B)', iinfo, n,
1096  $ jtype, ioldsd
1097  info = abs( iinfo )
1098  go to 210
1099  END IF
1100 *
1101  CALL dget52( .false., n, h, lda, t, lda, evectr, ldu,
1102  $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1103  result( 12 ) = dumma( 1 )
1104  IF( dumma( 2 ).GT.thresh ) THEN
1105  WRITE( nounit, fmt = 9998 )'Right', 'DTGEVC(HOWMNY=B)',
1106  $ dumma( 2 ), n, jtype, ioldsd
1107  END IF
1108 *
1109 * Tests 13--15 are done only on request
1110 *
1111  IF( tstdif ) THEN
1112 *
1113 * Do Tests 13--14
1114 *
1115  CALL dget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1116  $ work, result( 13 ) )
1117  CALL dget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1118  $ work, result( 14 ) )
1119 *
1120 * Do Test 15
1121 *
1122  temp1 = zero
1123  temp2 = zero
1124  DO 200 j = 1, n
1125  temp1 = max( temp1, abs( alphr1( j )-alphr3( j ) )+
1126  $ abs( alphi1( j )-alphi3( j ) ) )
1127  temp2 = max( temp2, abs( beta1( j )-beta3( j ) ) )
1128  200 continue
1129 *
1130  temp1 = temp1 / max( safmin, ulp*max( temp1, anorm ) )
1131  temp2 = temp2 / max( safmin, ulp*max( temp2, bnorm ) )
1132  result( 15 ) = max( temp1, temp2 )
1133  ntest = 15
1134  ELSE
1135  result( 13 ) = zero
1136  result( 14 ) = zero
1137  result( 15 ) = zero
1138  ntest = 12
1139  END IF
1140 *
1141 * End of Loop -- Check for RESULT(j) > THRESH
1142 *
1143  210 continue
1144 *
1145  ntestt = ntestt + ntest
1146 *
1147 * Print out tests which fail.
1148 *
1149  DO 220 jr = 1, ntest
1150  IF( result( jr ).GE.thresh ) THEN
1151 *
1152 * If this is the first test to fail,
1153 * print a header to the data file.
1154 *
1155  IF( nerrs.EQ.0 ) THEN
1156  WRITE( nounit, fmt = 9997 )'DGG'
1157 *
1158 * Matrix types
1159 *
1160  WRITE( nounit, fmt = 9996 )
1161  WRITE( nounit, fmt = 9995 )
1162  WRITE( nounit, fmt = 9994 )'Orthogonal'
1163 *
1164 * Tests performed
1165 *
1166  WRITE( nounit, fmt = 9993 )'orthogonal', '''',
1167  $ 'transpose', ( '''', j = 1, 10 )
1168 *
1169  END IF
1170  nerrs = nerrs + 1
1171  IF( result( jr ).LT.10000.0d0 ) THEN
1172  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1173  $ result( jr )
1174  ELSE
1175  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1176  $ result( jr )
1177  END IF
1178  END IF
1179  220 continue
1180 *
1181  230 continue
1182  240 continue
1183 *
1184 * Summary
1185 *
1186  CALL dlasum( 'DGG', nounit, nerrs, ntestt )
1187  return
1188 *
1189  9999 format( ' DCHKGG: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1190  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1191 *
1192  9998 format( ' DCHKGG: ', a, ' Eigenvectors from ', a, ' incorrectly ',
1193  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
1194  $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
1195  $ ')' )
1196 *
1197  9997 format( / 1x, a3, ' -- Real Generalized eigenvalue problem' )
1198 *
1199  9996 format( ' Matrix types (see DCHKGG for details): ' )
1200 *
1201  9995 format( ' Special Matrices:', 23x,
1202  $ '(J''=transposed Jordan block)',
1203  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1204  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
1205  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
1206  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
1207  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1208  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
1209  9994 format( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
1210  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
1211  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
1212  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
1213  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
1214  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
1215  $ '23=(small,large) 24=(small,small) 25=(large,large)',
1216  $ / ' 26=random O(1) matrices.' )
1217 *
1218  9993 format( / ' Tests performed: (H is Hessenberg, S is Schur, B, ',
1219  $ 'T, P are triangular,', / 20x, 'U, V, Q, and Z are ', a,
1220  $ ', l and r are the', / 20x,
1221  $ 'appropriate left and right eigenvectors, resp., a is',
1222  $ / 20x, 'alpha, b is beta, and ', a, ' means ', a, '.)',
1223  $ / ' 1 = | A - U H V', a,
1224  $ ' | / ( |A| n ulp ) 2 = | B - U T V', a,
1225  $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', a,
1226  $ ' | / ( n ulp ) 4 = | I - VV', a,
1227  $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', a,
1228  $ ' | / ( |H| n ulp )', 6x, '6 = | T - Q P Z', a,
1229  $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', a,
1230  $ ' | / ( n ulp ) 8 = | I - ZZ', a,
1231  $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', a,
1232  $ ' l | / const. 10 = max | ( b H - a T )', a,
1233  $ ' l | / const.', /
1234  $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1235  $ ' - a T ) r | / const.', / 1x )
1236 *
1237  9992 format( ' Matrix order=', i5, ', type=', i2, ', seed=',
1238  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
1239  9991 format( ' Matrix order=', i5, ', type=', i2, ', seed=',
1240  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
1241 *
1242 * End of DCHKGG
1243 *
1244  END