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