1 SUBROUTINE DDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
3 $ ALPHAI, BETA, WORK, LWORK, RESULT, BWORK,
4 $ INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * ), DOTYPE( * )
16 INTEGER ISEED( 4 ), NN( * )
17 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
18 $ B( LDA, * ), BETA( * ), Q( LDQ, * ),
19 $ RESULT( 13 ), S( LDA, * ), T( LDA, * ),
20 $ WORK( * ), Z( LDQ, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem driver DGGES.
28 *
29 * DGGES factors A and B as Q S Z' and Q T Z' , where ' means
30 * transpose, T is upper triangular, S is in generalized Schur form
31 * (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
32 * the 2x2 blocks corresponding to complex conjugate pairs of
33 * generalized eigenvalues), and Q and Z are orthogonal. It also
34 * computes the generalized eigenvalues (alpha(j),beta(j)), j=1,...,n,
35 * Thus, w(j) = alpha(j)/beta(j) is a root of the characteristic
36 * equation
37 * det( A - w(j) B ) = 0
38 * Optionally it also reorder the eigenvalues so that a selected
39 * cluster of eigenvalues appears in the leading diagonal block of the
40 * Schur forms.
41 *
42 * When DDRGES is called, a number of matrix "sizes" ("N's") and a
43 * number of matrix "TYPES" are specified. For each size ("N")
44 * and each TYPE of matrix, a pair of matrices (A, B) will be generated
45 * and used for testing. For each matrix pair, the following 13 tests
46 * will be performed and compared with the threshhold THRESH except
47 * the tests (5), (11) and (13).
48 *
49 *
50 * (1) | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
51 *
52 *
53 * (2) | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
54 *
55 *
56 * (3) | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
57 *
58 *
59 * (4) | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
60 *
61 * (5) if A is in Schur form (i.e. quasi-triangular form)
62 * (no sorting of eigenvalues)
63 *
64 * (6) if eigenvalues = diagonal blocks of the Schur form (S, T),
65 * i.e., test the maximum over j of D(j) where:
66 *
67 * if alpha(j) is real:
68 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
69 * D(j) = ------------------------ + -----------------------
70 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
71 *
72 * if alpha(j) is complex:
73 * | det( s S - w T ) |
74 * D(j) = ---------------------------------------------------
75 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
76 *
77 * and S and T are here the 2 x 2 diagonal blocks of S and T
78 * corresponding to the j-th and j+1-th eigenvalues.
79 * (no sorting of eigenvalues)
80 *
81 * (7) | (A,B) - Q (S,T) Z' | / ( | (A,B) | n ulp )
82 * (with sorting of eigenvalues).
83 *
84 * (8) | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
85 *
86 * (9) | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
87 *
88 * (10) if A is in Schur form (i.e. quasi-triangular form)
89 * (with sorting of eigenvalues).
90 *
91 * (11) if eigenvalues = diagonal blocks of the Schur form (S, T),
92 * i.e. test the maximum over j of D(j) where:
93 *
94 * if alpha(j) is real:
95 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
96 * D(j) = ------------------------ + -----------------------
97 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
98 *
99 * if alpha(j) is complex:
100 * | det( s S - w T ) |
101 * D(j) = ---------------------------------------------------
102 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
103 *
104 * and S and T are here the 2 x 2 diagonal blocks of S and T
105 * corresponding to the j-th and j+1-th eigenvalues.
106 * (with sorting of eigenvalues).
107 *
108 * (12) if sorting worked and SDIM is the number of eigenvalues
109 * which were SELECTed.
110 *
111 * Test Matrices
112 * =============
113 *
114 * The sizes of the test matrices are specified by an array
115 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
116 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
117 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
118 * Currently, the list of possible types is:
119 *
120 * (1) ( 0, 0 ) (a pair of zero matrices)
121 *
122 * (2) ( I, 0 ) (an identity and a zero matrix)
123 *
124 * (3) ( 0, I ) (an identity and a zero matrix)
125 *
126 * (4) ( I, I ) (a pair of identity matrices)
127 *
128 * t t
129 * (5) ( J , J ) (a pair of transposed Jordan blocks)
130 *
131 * t ( I 0 )
132 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
133 * ( 0 I ) ( 0 J )
134 * and I is a k x k identity and J a (k+1)x(k+1)
135 * Jordan block; k=(N-1)/2
136 *
137 * (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
138 * matrix with those diagonal entries.)
139 * (8) ( I, D )
140 *
141 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
142 *
143 * (10) ( small*D, big*I )
144 *
145 * (11) ( big*I, small*D )
146 *
147 * (12) ( small*I, big*D )
148 *
149 * (13) ( big*D, big*I )
150 *
151 * (14) ( small*D, small*I )
152 *
153 * (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
154 * D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
155 * t t
156 * (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
157 *
158 * (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
159 * with random O(1) entries above the diagonal
160 * and diagonal entries diag(T1) =
161 * ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
162 * ( 0, N-3, N-4,..., 1, 0, 0 )
163 *
164 * (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
165 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
166 * s = machine precision.
167 *
168 * (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
169 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
170 *
171 * N-5
172 * (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
173 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
174 *
175 * (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
176 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
177 * where r1,..., r(N-4) are random.
178 *
179 * (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
180 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
181 *
182 * (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
183 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
184 *
185 * (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
186 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
187 *
188 * (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
189 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
190 *
191 * (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
192 * matrices.
193 *
194 *
195 * Arguments
196 * =========
197 *
198 * NSIZES (input) INTEGER
199 * The number of sizes of matrices to use. If it is zero,
200 * DDRGES does nothing. NSIZES >= 0.
201 *
202 * NN (input) INTEGER array, dimension (NSIZES)
203 * An array containing the sizes to be used for the matrices.
204 * Zero values will be skipped. NN >= 0.
205 *
206 * NTYPES (input) INTEGER
207 * The number of elements in DOTYPE. If it is zero, DDRGES
208 * does nothing. It must be at least zero. If it is MAXTYP+1
209 * and NSIZES is 1, then an additional type, MAXTYP+1 is
210 * defined, which is to use whatever matrix is in A on input.
211 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
212 * DOTYPE(MAXTYP+1) is .TRUE. .
213 *
214 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
215 * If DOTYPE(j) is .TRUE., then for each size in NN a
216 * matrix of that size and of type j will be generated.
217 * If NTYPES is smaller than the maximum number of types
218 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
219 * MAXTYP will not be generated. If NTYPES is larger
220 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
221 * will be ignored.
222 *
223 * ISEED (input/output) INTEGER array, dimension (4)
224 * On entry ISEED specifies the seed of the random number
225 * generator. The array elements should be between 0 and 4095;
226 * if not they will be reduced mod 4096. Also, ISEED(4) must
227 * be odd. The random number generator uses a linear
228 * congruential sequence limited to small integers, and so
229 * should produce machine independent random numbers. The
230 * values of ISEED are changed on exit, and can be used in the
231 * next call to DDRGES to continue the same random number
232 * sequence.
233 *
234 * THRESH (input) DOUBLE PRECISION
235 * A test will count as "failed" if the "error", computed as
236 * described above, exceeds THRESH. Note that the error is
237 * scaled to be O(1), so THRESH should be a reasonably small
238 * multiple of 1, e.g., 10 or 100. In particular, it should
239 * not depend on the precision (single vs. double) or the size
240 * of the matrix. THRESH >= 0.
241 *
242 * NOUNIT (input) INTEGER
243 * The FORTRAN unit number for printing out error messages
244 * (e.g., if a routine returns IINFO not equal to 0.)
245 *
246 * A (input/workspace) DOUBLE PRECISION array,
247 * dimension(LDA, max(NN))
248 * Used to hold the original A matrix. Used as input only
249 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
250 * DOTYPE(MAXTYP+1)=.TRUE.
251 *
252 * LDA (input) INTEGER
253 * The leading dimension of A, B, S, and T.
254 * It must be at least 1 and at least max( NN ).
255 *
256 * B (input/workspace) DOUBLE PRECISION array,
257 * dimension(LDA, max(NN))
258 * Used to hold the original B matrix. Used as input only
259 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
260 * DOTYPE(MAXTYP+1)=.TRUE.
261 *
262 * S (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
263 * The Schur form matrix computed from A by DGGES. On exit, S
264 * contains the Schur form matrix corresponding to the matrix
265 * in A.
266 *
267 * T (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
268 * The upper triangular matrix computed from B by DGGES.
269 *
270 * Q (workspace) DOUBLE PRECISION array, dimension (LDQ, max(NN))
271 * The (left) orthogonal matrix computed by DGGES.
272 *
273 * LDQ (input) INTEGER
274 * The leading dimension of Q and Z. It must
275 * be at least 1 and at least max( NN ).
276 *
277 * Z (workspace) DOUBLE PRECISION array, dimension( LDQ, max(NN) )
278 * The (right) orthogonal matrix computed by DGGES.
279 *
280 * ALPHAR (workspace) DOUBLE PRECISION array, dimension (max(NN))
281 * ALPHAI (workspace) DOUBLE PRECISION array, dimension (max(NN))
282 * BETA (workspace) DOUBLE PRECISION array, dimension (max(NN))
283 * The generalized eigenvalues of (A,B) computed by DGGES.
284 * ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
285 * generalized eigenvalue of A and B.
286 *
287 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
288 *
289 * LWORK (input) INTEGER
290 * The dimension of the array WORK.
291 * LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
292 * matrix dimension.
293 *
294 * RESULT (output) DOUBLE PRECISION array, dimension (15)
295 * The values computed by the tests described above.
296 * The values are currently limited to 1/ulp, to avoid overflow.
297 *
298 * BWORK (workspace) LOGICAL array, dimension (N)
299 *
300 * INFO (output) INTEGER
301 * = 0: successful exit
302 * < 0: if INFO = -i, the i-th argument had an illegal value.
303 * > 0: A routine returned an error code. INFO is the
304 * absolute value of the INFO value returned.
305 *
306 * =====================================================================
307 *
308 * .. Parameters ..
309 DOUBLE PRECISION ZERO, ONE
310 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
311 INTEGER MAXTYP
312 PARAMETER ( MAXTYP = 26 )
313 * ..
314 * .. Local Scalars ..
315 LOGICAL BADNN, ILABAD
316 CHARACTER SORT
317 INTEGER I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
318 $ JSIZE, JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES,
319 $ N, N1, NB, NERRS, NMATS, NMAX, NTEST, NTESTT,
320 $ RSUB, SDIM
321 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
322 * ..
323 * .. Local Arrays ..
324 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
325 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
326 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
327 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
328 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
329 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
330 DOUBLE PRECISION RMAGN( 0: 3 )
331 * ..
332 * .. External Functions ..
333 LOGICAL DLCTES
334 INTEGER ILAENV
335 DOUBLE PRECISION DLAMCH, DLARND
336 EXTERNAL DLCTES, ILAENV, DLAMCH, DLARND
337 * ..
338 * .. External Subroutines ..
339 EXTERNAL ALASVM, DGET51, DGET53, DGET54, DGGES, DLABAD,
340 $ DLACPY, DLARFG, DLASET, DLATM4, DORM2R, XERBLA
341 * ..
342 * .. Intrinsic Functions ..
343 INTRINSIC ABS, DBLE, MAX, MIN, SIGN
344 * ..
345 * .. Data statements ..
346 DATA KCLASS / 15*1, 10*2, 1*3 /
347 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
348 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
349 DATA KADD / 0, 0, 0, 0, 3, 2 /
350 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
351 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
352 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
353 $ 1, 1, -4, 2, -4, 8*8, 0 /
354 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
355 $ 4*5, 4*3, 1 /
356 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
357 $ 4*6, 4*4, 1 /
358 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
359 $ 2, 1 /
360 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
361 $ 2, 1 /
362 DATA KTRIAN / 16*0, 10*1 /
363 DATA IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
364 $ 5*2, 0 /
365 DATA IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
366 * ..
367 * .. Executable Statements ..
368 *
369 * Check for errors
370 *
371 INFO = 0
372 *
373 BADNN = .FALSE.
374 NMAX = 1
375 DO 10 J = 1, NSIZES
376 NMAX = MAX( NMAX, NN( J ) )
377 IF( NN( J ).LT.0 )
378 $ BADNN = .TRUE.
379 10 CONTINUE
380 *
381 IF( NSIZES.LT.0 ) THEN
382 INFO = -1
383 ELSE IF( BADNN ) THEN
384 INFO = -2
385 ELSE IF( NTYPES.LT.0 ) THEN
386 INFO = -3
387 ELSE IF( THRESH.LT.ZERO ) THEN
388 INFO = -6
389 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
390 INFO = -9
391 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
392 INFO = -14
393 END IF
394 *
395 * Compute workspace
396 * (Note: Comments in the code beginning "Workspace:" describe the
397 * minimal amount of workspace needed at that point in the code,
398 * as well as the preferred amount for good performance.
399 * NB refers to the optimal block size for the immediately
400 * following subroutine, as returned by ILAENV.
401 *
402 MINWRK = 1
403 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
404 MINWRK = MAX( 10*( NMAX+1 ), 3*NMAX*NMAX )
405 NB = MAX( 1, ILAENV( 1, 'DGEQRF', ' ', NMAX, NMAX, -1, -1 ),
406 $ ILAENV( 1, 'DORMQR', 'LT', NMAX, NMAX, NMAX, -1 ),
407 $ ILAENV( 1, 'DORGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
408 MAXWRK = MAX( 10*( NMAX+1 ), 2*NMAX+NMAX*NB, 3*NMAX*NMAX )
409 WORK( 1 ) = MAXWRK
410 END IF
411 *
412 IF( LWORK.LT.MINWRK )
413 $ INFO = -20
414 *
415 IF( INFO.NE.0 ) THEN
416 CALL XERBLA( 'DDRGES', -INFO )
417 RETURN
418 END IF
419 *
420 * Quick return if possible
421 *
422 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
423 $ RETURN
424 *
425 SAFMIN = DLAMCH( 'Safe minimum' )
426 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
427 SAFMIN = SAFMIN / ULP
428 SAFMAX = ONE / SAFMIN
429 CALL DLABAD( SAFMIN, SAFMAX )
430 ULPINV = ONE / ULP
431 *
432 * The values RMAGN(2:3) depend on N, see below.
433 *
434 RMAGN( 0 ) = ZERO
435 RMAGN( 1 ) = ONE
436 *
437 * Loop over matrix sizes
438 *
439 NTESTT = 0
440 NERRS = 0
441 NMATS = 0
442 *
443 DO 190 JSIZE = 1, NSIZES
444 N = NN( JSIZE )
445 N1 = MAX( 1, N )
446 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
447 RMAGN( 3 ) = SAFMIN*ULPINV*DBLE( N1 )
448 *
449 IF( NSIZES.NE.1 ) THEN
450 MTYPES = MIN( MAXTYP, NTYPES )
451 ELSE
452 MTYPES = MIN( MAXTYP+1, NTYPES )
453 END IF
454 *
455 * Loop over matrix types
456 *
457 DO 180 JTYPE = 1, MTYPES
458 IF( .NOT.DOTYPE( JTYPE ) )
459 $ GO TO 180
460 NMATS = NMATS + 1
461 NTEST = 0
462 *
463 * Save ISEED in case of an error.
464 *
465 DO 20 J = 1, 4
466 IOLDSD( J ) = ISEED( J )
467 20 CONTINUE
468 *
469 * Initialize RESULT
470 *
471 DO 30 J = 1, 13
472 RESULT( J ) = ZERO
473 30 CONTINUE
474 *
475 * Generate test matrices A and B
476 *
477 * Description of control parameters:
478 *
479 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
480 * =3 means random.
481 * KATYPE: the "type" to be passed to DLATM4 for computing A.
482 * KAZERO: the pattern of zeros on the diagonal for A:
483 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
484 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
485 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
486 * non-zero entries.)
487 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
488 * =2: large, =3: small.
489 * IASIGN: 1 if the diagonal elements of A are to be
490 * multiplied by a random magnitude 1 number, =2 if
491 * randomly chosen diagonal blocks are to be rotated
492 * to form 2x2 blocks.
493 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
494 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
495 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
496 * RMAGN: used to implement KAMAGN and KBMAGN.
497 *
498 IF( MTYPES.GT.MAXTYP )
499 $ GO TO 110
500 IINFO = 0
501 IF( KCLASS( JTYPE ).LT.3 ) THEN
502 *
503 * Generate A (w/o rotation)
504 *
505 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
506 IN = 2*( ( N-1 ) / 2 ) + 1
507 IF( IN.NE.N )
508 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
509 ELSE
510 IN = N
511 END IF
512 CALL DLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
513 $ KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ),
514 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
515 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
516 $ ISEED, A, LDA )
517 IADD = KADD( KAZERO( JTYPE ) )
518 IF( IADD.GT.0 .AND. IADD.LE.N )
519 $ A( IADD, IADD ) = ONE
520 *
521 * Generate B (w/o rotation)
522 *
523 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
524 IN = 2*( ( N-1 ) / 2 ) + 1
525 IF( IN.NE.N )
526 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDA )
527 ELSE
528 IN = N
529 END IF
530 CALL DLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
531 $ KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ),
532 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
533 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
534 $ ISEED, B, LDA )
535 IADD = KADD( KBZERO( JTYPE ) )
536 IF( IADD.NE.0 .AND. IADD.LE.N )
537 $ B( IADD, IADD ) = ONE
538 *
539 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
540 *
541 * Include rotations
542 *
543 * Generate Q, Z as Householder transformations times
544 * a diagonal matrix.
545 *
546 DO 50 JC = 1, N - 1
547 DO 40 JR = JC, N
548 Q( JR, JC ) = DLARND( 3, ISEED )
549 Z( JR, JC ) = DLARND( 3, ISEED )
550 40 CONTINUE
551 CALL DLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
552 $ WORK( JC ) )
553 WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) )
554 Q( JC, JC ) = ONE
555 CALL DLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
556 $ WORK( N+JC ) )
557 WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) )
558 Z( JC, JC ) = ONE
559 50 CONTINUE
560 Q( N, N ) = ONE
561 WORK( N ) = ZERO
562 WORK( 3*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
563 Z( N, N ) = ONE
564 WORK( 2*N ) = ZERO
565 WORK( 4*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
566 *
567 * Apply the diagonal matrices
568 *
569 DO 70 JC = 1, N
570 DO 60 JR = 1, N
571 A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
572 $ A( JR, JC )
573 B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
574 $ B( JR, JC )
575 60 CONTINUE
576 70 CONTINUE
577 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
578 $ LDA, WORK( 2*N+1 ), IINFO )
579 IF( IINFO.NE.0 )
580 $ GO TO 100
581 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
582 $ A, LDA, WORK( 2*N+1 ), IINFO )
583 IF( IINFO.NE.0 )
584 $ GO TO 100
585 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
586 $ LDA, WORK( 2*N+1 ), IINFO )
587 IF( IINFO.NE.0 )
588 $ GO TO 100
589 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
590 $ B, LDA, WORK( 2*N+1 ), IINFO )
591 IF( IINFO.NE.0 )
592 $ GO TO 100
593 END IF
594 ELSE
595 *
596 * Random matrices
597 *
598 DO 90 JC = 1, N
599 DO 80 JR = 1, N
600 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
601 $ DLARND( 2, ISEED )
602 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
603 $ DLARND( 2, ISEED )
604 80 CONTINUE
605 90 CONTINUE
606 END IF
607 *
608 100 CONTINUE
609 *
610 IF( IINFO.NE.0 ) THEN
611 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
612 $ IOLDSD
613 INFO = ABS( IINFO )
614 RETURN
615 END IF
616 *
617 110 CONTINUE
618 *
619 DO 120 I = 1, 13
620 RESULT( I ) = -ONE
621 120 CONTINUE
622 *
623 * Test with and without sorting of eigenvalues
624 *
625 DO 150 ISORT = 0, 1
626 IF( ISORT.EQ.0 ) THEN
627 SORT = 'N'
628 RSUB = 0
629 ELSE
630 SORT = 'S'
631 RSUB = 5
632 END IF
633 *
634 * Call DGGES to compute H, T, Q, Z, alpha, and beta.
635 *
636 CALL DLACPY( 'Full', N, N, A, LDA, S, LDA )
637 CALL DLACPY( 'Full', N, N, B, LDA, T, LDA )
638 NTEST = 1 + RSUB + ISORT
639 RESULT( 1+RSUB+ISORT ) = ULPINV
640 CALL DGGES( 'V', 'V', SORT, DLCTES, N, S, LDA, T, LDA,
641 $ SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDQ,
642 $ WORK, LWORK, BWORK, IINFO )
643 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
644 RESULT( 1+RSUB+ISORT ) = ULPINV
645 WRITE( NOUNIT, FMT = 9999 )'DGGES', IINFO, N, JTYPE,
646 $ IOLDSD
647 INFO = ABS( IINFO )
648 GO TO 160
649 END IF
650 *
651 NTEST = 4 + RSUB
652 *
653 * Do tests 1--4 (or tests 7--9 when reordering )
654 *
655 IF( ISORT.EQ.0 ) THEN
656 CALL DGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ,
657 $ WORK, RESULT( 1 ) )
658 CALL DGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ,
659 $ WORK, RESULT( 2 ) )
660 ELSE
661 CALL DGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q,
662 $ LDQ, Z, LDQ, WORK, RESULT( 7 ) )
663 END IF
664 CALL DGET51( 3, N, A, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
665 $ RESULT( 3+RSUB ) )
666 CALL DGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
667 $ RESULT( 4+RSUB ) )
668 *
669 * Do test 5 and 6 (or Tests 10 and 11 when reordering):
670 * check Schur form of A and compare eigenvalues with
671 * diagonals.
672 *
673 NTEST = 6 + RSUB
674 TEMP1 = ZERO
675 *
676 DO 130 J = 1, N
677 ILABAD = .FALSE.
678 IF( ALPHAI( J ).EQ.ZERO ) THEN
679 TEMP2 = ( ABS( ALPHAR( J )-S( J, J ) ) /
680 $ MAX( SAFMIN, ABS( ALPHAR( J ) ), ABS( S( J,
681 $ J ) ) )+ABS( BETA( J )-T( J, J ) ) /
682 $ MAX( SAFMIN, ABS( BETA( J ) ), ABS( T( J,
683 $ J ) ) ) ) / ULP
684 *
685 IF( J.LT.N ) THEN
686 IF( S( J+1, J ).NE.ZERO ) THEN
687 ILABAD = .TRUE.
688 RESULT( 5+RSUB ) = ULPINV
689 END IF
690 END IF
691 IF( J.GT.1 ) THEN
692 IF( S( J, J-1 ).NE.ZERO ) THEN
693 ILABAD = .TRUE.
694 RESULT( 5+RSUB ) = ULPINV
695 END IF
696 END IF
697 *
698 ELSE
699 IF( ALPHAI( J ).GT.ZERO ) THEN
700 I1 = J
701 ELSE
702 I1 = J - 1
703 END IF
704 IF( I1.LE.0 .OR. I1.GE.N ) THEN
705 ILABAD = .TRUE.
706 ELSE IF( I1.LT.N-1 ) THEN
707 IF( S( I1+2, I1+1 ).NE.ZERO ) THEN
708 ILABAD = .TRUE.
709 RESULT( 5+RSUB ) = ULPINV
710 END IF
711 ELSE IF( I1.GT.1 ) THEN
712 IF( S( I1, I1-1 ).NE.ZERO ) THEN
713 ILABAD = .TRUE.
714 RESULT( 5+RSUB ) = ULPINV
715 END IF
716 END IF
717 IF( .NOT.ILABAD ) THEN
718 CALL DGET53( S( I1, I1 ), LDA, T( I1, I1 ), LDA,
719 $ BETA( J ), ALPHAR( J ),
720 $ ALPHAI( J ), TEMP2, IERR )
721 IF( IERR.GE.3 ) THEN
722 WRITE( NOUNIT, FMT = 9998 )IERR, J, N,
723 $ JTYPE, IOLDSD
724 INFO = ABS( IERR )
725 END IF
726 ELSE
727 TEMP2 = ULPINV
728 END IF
729 *
730 END IF
731 TEMP1 = MAX( TEMP1, TEMP2 )
732 IF( ILABAD ) THEN
733 WRITE( NOUNIT, FMT = 9997 )J, N, JTYPE, IOLDSD
734 END IF
735 130 CONTINUE
736 RESULT( 6+RSUB ) = TEMP1
737 *
738 IF( ISORT.GE.1 ) THEN
739 *
740 * Do test 12
741 *
742 NTEST = 12
743 RESULT( 12 ) = ZERO
744 KNTEIG = 0
745 DO 140 I = 1, N
746 IF( DLCTES( ALPHAR( I ), ALPHAI( I ),
747 $ BETA( I ) ) .OR. DLCTES( ALPHAR( I ),
748 $ -ALPHAI( I ), BETA( I ) ) ) THEN
749 KNTEIG = KNTEIG + 1
750 END IF
751 IF( I.LT.N ) THEN
752 IF( ( DLCTES( ALPHAR( I+1 ), ALPHAI( I+1 ),
753 $ BETA( I+1 ) ) .OR. DLCTES( ALPHAR( I+1 ),
754 $ -ALPHAI( I+1 ), BETA( I+1 ) ) ) .AND.
755 $ ( .NOT.( DLCTES( ALPHAR( I ), ALPHAI( I ),
756 $ BETA( I ) ) .OR. DLCTES( ALPHAR( I ),
757 $ -ALPHAI( I ), BETA( I ) ) ) ) .AND.
758 $ IINFO.NE.N+2 ) THEN
759 RESULT( 12 ) = ULPINV
760 END IF
761 END IF
762 140 CONTINUE
763 IF( SDIM.NE.KNTEIG ) THEN
764 RESULT( 12 ) = ULPINV
765 END IF
766 END IF
767 *
768 150 CONTINUE
769 *
770 * End of Loop -- Check for RESULT(j) > THRESH
771 *
772 160 CONTINUE
773 *
774 NTESTT = NTESTT + NTEST
775 *
776 * Print out tests which fail.
777 *
778 DO 170 JR = 1, NTEST
779 IF( RESULT( JR ).GE.THRESH ) THEN
780 *
781 * If this is the first test to fail,
782 * print a header to the data file.
783 *
784 IF( NERRS.EQ.0 ) THEN
785 WRITE( NOUNIT, FMT = 9996 )'DGS'
786 *
787 * Matrix types
788 *
789 WRITE( NOUNIT, FMT = 9995 )
790 WRITE( NOUNIT, FMT = 9994 )
791 WRITE( NOUNIT, FMT = 9993 )'Orthogonal'
792 *
793 * Tests performed
794 *
795 WRITE( NOUNIT, FMT = 9992 )'orthogonal', '''',
796 $ 'transpose', ( '''', J = 1, 8 )
797 *
798 END IF
799 NERRS = NERRS + 1
800 IF( RESULT( JR ).LT.10000.0D0 ) THEN
801 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
802 $ RESULT( JR )
803 ELSE
804 WRITE( NOUNIT, FMT = 9990 )N, JTYPE, IOLDSD, JR,
805 $ RESULT( JR )
806 END IF
807 END IF
808 170 CONTINUE
809 *
810 180 CONTINUE
811 190 CONTINUE
812 *
813 * Summary
814 *
815 CALL ALASVM( 'DGS', NOUNIT, NERRS, NTESTT, 0 )
816 *
817 WORK( 1 ) = MAXWRK
818 *
819 RETURN
820 *
821 9999 FORMAT( ' DDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
822 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
823 *
824 9998 FORMAT( ' DDRGES: DGET53 returned INFO=', I1, ' for eigenvalue ',
825 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(',
826 $ 4( I4, ',' ), I5, ')' )
827 *
828 9997 FORMAT( ' DDRGES: S not in Schur form at eigenvalue ', I6, '.',
829 $ / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
830 $ I5, ')' )
831 *
832 9996 FORMAT( / 1X, A3, ' -- Real Generalized Schur form driver' )
833 *
834 9995 FORMAT( ' Matrix types (see DDRGES for details): ' )
835 *
836 9994 FORMAT( ' Special Matrices:', 23X,
837 $ '(J''=transposed Jordan block)',
838 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
839 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
840 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
841 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
842 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
843 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
844 9993 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
845 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
846 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
847 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
848 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
849 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
850 $ '23=(small,large) 24=(small,small) 25=(large,large)',
851 $ / ' 26=random O(1) matrices.' )
852 *
853 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
854 $ 'Q and Z are ', A, ',', / 19X,
855 $ 'l and r are the appropriate left and right', / 19X,
856 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
857 $ ' means ', A, '.)', / ' Without ordering: ',
858 $ / ' 1 = | A - Q S Z', A,
859 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
860 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
861 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
862 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
863 $ / ' 6 = difference between (alpha,beta)',
864 $ ' and diagonals of (S,T)', / ' With ordering: ',
865 $ / ' 7 = | (A,B) - Q (S,T) Z', A,
866 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', A,
867 $ ' | / ( n ulp ) 9 = | I - ZZ', A,
868 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
869 $ / ' 11 = difference between (alpha,beta) and diagonals',
870 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
871 $ 'selected eigenvalues', / )
872 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
873 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
874 9990 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
875 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
876 *
877 * End of DDRGES
878 *
879 END
2 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
3 $ ALPHAI, BETA, WORK, LWORK, RESULT, BWORK,
4 $ INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
12 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * ), DOTYPE( * )
16 INTEGER ISEED( 4 ), NN( * )
17 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
18 $ B( LDA, * ), BETA( * ), Q( LDQ, * ),
19 $ RESULT( 13 ), S( LDA, * ), T( LDA, * ),
20 $ WORK( * ), Z( LDQ, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem driver DGGES.
28 *
29 * DGGES factors A and B as Q S Z' and Q T Z' , where ' means
30 * transpose, T is upper triangular, S is in generalized Schur form
31 * (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
32 * the 2x2 blocks corresponding to complex conjugate pairs of
33 * generalized eigenvalues), and Q and Z are orthogonal. It also
34 * computes the generalized eigenvalues (alpha(j),beta(j)), j=1,...,n,
35 * Thus, w(j) = alpha(j)/beta(j) is a root of the characteristic
36 * equation
37 * det( A - w(j) B ) = 0
38 * Optionally it also reorder the eigenvalues so that a selected
39 * cluster of eigenvalues appears in the leading diagonal block of the
40 * Schur forms.
41 *
42 * When DDRGES is called, a number of matrix "sizes" ("N's") and a
43 * number of matrix "TYPES" are specified. For each size ("N")
44 * and each TYPE of matrix, a pair of matrices (A, B) will be generated
45 * and used for testing. For each matrix pair, the following 13 tests
46 * will be performed and compared with the threshhold THRESH except
47 * the tests (5), (11) and (13).
48 *
49 *
50 * (1) | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
51 *
52 *
53 * (2) | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
54 *
55 *
56 * (3) | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
57 *
58 *
59 * (4) | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
60 *
61 * (5) if A is in Schur form (i.e. quasi-triangular form)
62 * (no sorting of eigenvalues)
63 *
64 * (6) if eigenvalues = diagonal blocks of the Schur form (S, T),
65 * i.e., test the maximum over j of D(j) where:
66 *
67 * if alpha(j) is real:
68 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
69 * D(j) = ------------------------ + -----------------------
70 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
71 *
72 * if alpha(j) is complex:
73 * | det( s S - w T ) |
74 * D(j) = ---------------------------------------------------
75 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
76 *
77 * and S and T are here the 2 x 2 diagonal blocks of S and T
78 * corresponding to the j-th and j+1-th eigenvalues.
79 * (no sorting of eigenvalues)
80 *
81 * (7) | (A,B) - Q (S,T) Z' | / ( | (A,B) | n ulp )
82 * (with sorting of eigenvalues).
83 *
84 * (8) | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
85 *
86 * (9) | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
87 *
88 * (10) if A is in Schur form (i.e. quasi-triangular form)
89 * (with sorting of eigenvalues).
90 *
91 * (11) if eigenvalues = diagonal blocks of the Schur form (S, T),
92 * i.e. test the maximum over j of D(j) where:
93 *
94 * if alpha(j) is real:
95 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
96 * D(j) = ------------------------ + -----------------------
97 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
98 *
99 * if alpha(j) is complex:
100 * | det( s S - w T ) |
101 * D(j) = ---------------------------------------------------
102 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
103 *
104 * and S and T are here the 2 x 2 diagonal blocks of S and T
105 * corresponding to the j-th and j+1-th eigenvalues.
106 * (with sorting of eigenvalues).
107 *
108 * (12) if sorting worked and SDIM is the number of eigenvalues
109 * which were SELECTed.
110 *
111 * Test Matrices
112 * =============
113 *
114 * The sizes of the test matrices are specified by an array
115 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
116 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
117 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
118 * Currently, the list of possible types is:
119 *
120 * (1) ( 0, 0 ) (a pair of zero matrices)
121 *
122 * (2) ( I, 0 ) (an identity and a zero matrix)
123 *
124 * (3) ( 0, I ) (an identity and a zero matrix)
125 *
126 * (4) ( I, I ) (a pair of identity matrices)
127 *
128 * t t
129 * (5) ( J , J ) (a pair of transposed Jordan blocks)
130 *
131 * t ( I 0 )
132 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
133 * ( 0 I ) ( 0 J )
134 * and I is a k x k identity and J a (k+1)x(k+1)
135 * Jordan block; k=(N-1)/2
136 *
137 * (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
138 * matrix with those diagonal entries.)
139 * (8) ( I, D )
140 *
141 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
142 *
143 * (10) ( small*D, big*I )
144 *
145 * (11) ( big*I, small*D )
146 *
147 * (12) ( small*I, big*D )
148 *
149 * (13) ( big*D, big*I )
150 *
151 * (14) ( small*D, small*I )
152 *
153 * (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
154 * D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
155 * t t
156 * (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
157 *
158 * (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
159 * with random O(1) entries above the diagonal
160 * and diagonal entries diag(T1) =
161 * ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
162 * ( 0, N-3, N-4,..., 1, 0, 0 )
163 *
164 * (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
165 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
166 * s = machine precision.
167 *
168 * (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
169 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
170 *
171 * N-5
172 * (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
173 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
174 *
175 * (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
176 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
177 * where r1,..., r(N-4) are random.
178 *
179 * (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
180 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
181 *
182 * (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
183 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
184 *
185 * (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
186 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
187 *
188 * (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
189 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
190 *
191 * (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
192 * matrices.
193 *
194 *
195 * Arguments
196 * =========
197 *
198 * NSIZES (input) INTEGER
199 * The number of sizes of matrices to use. If it is zero,
200 * DDRGES does nothing. NSIZES >= 0.
201 *
202 * NN (input) INTEGER array, dimension (NSIZES)
203 * An array containing the sizes to be used for the matrices.
204 * Zero values will be skipped. NN >= 0.
205 *
206 * NTYPES (input) INTEGER
207 * The number of elements in DOTYPE. If it is zero, DDRGES
208 * does nothing. It must be at least zero. If it is MAXTYP+1
209 * and NSIZES is 1, then an additional type, MAXTYP+1 is
210 * defined, which is to use whatever matrix is in A on input.
211 * This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
212 * DOTYPE(MAXTYP+1) is .TRUE. .
213 *
214 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
215 * If DOTYPE(j) is .TRUE., then for each size in NN a
216 * matrix of that size and of type j will be generated.
217 * If NTYPES is smaller than the maximum number of types
218 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
219 * MAXTYP will not be generated. If NTYPES is larger
220 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
221 * will be ignored.
222 *
223 * ISEED (input/output) INTEGER array, dimension (4)
224 * On entry ISEED specifies the seed of the random number
225 * generator. The array elements should be between 0 and 4095;
226 * if not they will be reduced mod 4096. Also, ISEED(4) must
227 * be odd. The random number generator uses a linear
228 * congruential sequence limited to small integers, and so
229 * should produce machine independent random numbers. The
230 * values of ISEED are changed on exit, and can be used in the
231 * next call to DDRGES to continue the same random number
232 * sequence.
233 *
234 * THRESH (input) DOUBLE PRECISION
235 * A test will count as "failed" if the "error", computed as
236 * described above, exceeds THRESH. Note that the error is
237 * scaled to be O(1), so THRESH should be a reasonably small
238 * multiple of 1, e.g., 10 or 100. In particular, it should
239 * not depend on the precision (single vs. double) or the size
240 * of the matrix. THRESH >= 0.
241 *
242 * NOUNIT (input) INTEGER
243 * The FORTRAN unit number for printing out error messages
244 * (e.g., if a routine returns IINFO not equal to 0.)
245 *
246 * A (input/workspace) DOUBLE PRECISION array,
247 * dimension(LDA, max(NN))
248 * Used to hold the original A matrix. Used as input only
249 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
250 * DOTYPE(MAXTYP+1)=.TRUE.
251 *
252 * LDA (input) INTEGER
253 * The leading dimension of A, B, S, and T.
254 * It must be at least 1 and at least max( NN ).
255 *
256 * B (input/workspace) DOUBLE PRECISION array,
257 * dimension(LDA, max(NN))
258 * Used to hold the original B matrix. Used as input only
259 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
260 * DOTYPE(MAXTYP+1)=.TRUE.
261 *
262 * S (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
263 * The Schur form matrix computed from A by DGGES. On exit, S
264 * contains the Schur form matrix corresponding to the matrix
265 * in A.
266 *
267 * T (workspace) DOUBLE PRECISION array, dimension (LDA, max(NN))
268 * The upper triangular matrix computed from B by DGGES.
269 *
270 * Q (workspace) DOUBLE PRECISION array, dimension (LDQ, max(NN))
271 * The (left) orthogonal matrix computed by DGGES.
272 *
273 * LDQ (input) INTEGER
274 * The leading dimension of Q and Z. It must
275 * be at least 1 and at least max( NN ).
276 *
277 * Z (workspace) DOUBLE PRECISION array, dimension( LDQ, max(NN) )
278 * The (right) orthogonal matrix computed by DGGES.
279 *
280 * ALPHAR (workspace) DOUBLE PRECISION array, dimension (max(NN))
281 * ALPHAI (workspace) DOUBLE PRECISION array, dimension (max(NN))
282 * BETA (workspace) DOUBLE PRECISION array, dimension (max(NN))
283 * The generalized eigenvalues of (A,B) computed by DGGES.
284 * ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
285 * generalized eigenvalue of A and B.
286 *
287 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
288 *
289 * LWORK (input) INTEGER
290 * The dimension of the array WORK.
291 * LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
292 * matrix dimension.
293 *
294 * RESULT (output) DOUBLE PRECISION array, dimension (15)
295 * The values computed by the tests described above.
296 * The values are currently limited to 1/ulp, to avoid overflow.
297 *
298 * BWORK (workspace) LOGICAL array, dimension (N)
299 *
300 * INFO (output) INTEGER
301 * = 0: successful exit
302 * < 0: if INFO = -i, the i-th argument had an illegal value.
303 * > 0: A routine returned an error code. INFO is the
304 * absolute value of the INFO value returned.
305 *
306 * =====================================================================
307 *
308 * .. Parameters ..
309 DOUBLE PRECISION ZERO, ONE
310 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
311 INTEGER MAXTYP
312 PARAMETER ( MAXTYP = 26 )
313 * ..
314 * .. Local Scalars ..
315 LOGICAL BADNN, ILABAD
316 CHARACTER SORT
317 INTEGER I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
318 $ JSIZE, JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES,
319 $ N, N1, NB, NERRS, NMATS, NMAX, NTEST, NTESTT,
320 $ RSUB, SDIM
321 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
322 * ..
323 * .. Local Arrays ..
324 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
325 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
326 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
327 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
328 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
329 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
330 DOUBLE PRECISION RMAGN( 0: 3 )
331 * ..
332 * .. External Functions ..
333 LOGICAL DLCTES
334 INTEGER ILAENV
335 DOUBLE PRECISION DLAMCH, DLARND
336 EXTERNAL DLCTES, ILAENV, DLAMCH, DLARND
337 * ..
338 * .. External Subroutines ..
339 EXTERNAL ALASVM, DGET51, DGET53, DGET54, DGGES, DLABAD,
340 $ DLACPY, DLARFG, DLASET, DLATM4, DORM2R, XERBLA
341 * ..
342 * .. Intrinsic Functions ..
343 INTRINSIC ABS, DBLE, MAX, MIN, SIGN
344 * ..
345 * .. Data statements ..
346 DATA KCLASS / 15*1, 10*2, 1*3 /
347 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
348 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
349 DATA KADD / 0, 0, 0, 0, 3, 2 /
350 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
351 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
352 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
353 $ 1, 1, -4, 2, -4, 8*8, 0 /
354 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
355 $ 4*5, 4*3, 1 /
356 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
357 $ 4*6, 4*4, 1 /
358 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
359 $ 2, 1 /
360 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
361 $ 2, 1 /
362 DATA KTRIAN / 16*0, 10*1 /
363 DATA IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
364 $ 5*2, 0 /
365 DATA IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
366 * ..
367 * .. Executable Statements ..
368 *
369 * Check for errors
370 *
371 INFO = 0
372 *
373 BADNN = .FALSE.
374 NMAX = 1
375 DO 10 J = 1, NSIZES
376 NMAX = MAX( NMAX, NN( J ) )
377 IF( NN( J ).LT.0 )
378 $ BADNN = .TRUE.
379 10 CONTINUE
380 *
381 IF( NSIZES.LT.0 ) THEN
382 INFO = -1
383 ELSE IF( BADNN ) THEN
384 INFO = -2
385 ELSE IF( NTYPES.LT.0 ) THEN
386 INFO = -3
387 ELSE IF( THRESH.LT.ZERO ) THEN
388 INFO = -6
389 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
390 INFO = -9
391 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
392 INFO = -14
393 END IF
394 *
395 * Compute workspace
396 * (Note: Comments in the code beginning "Workspace:" describe the
397 * minimal amount of workspace needed at that point in the code,
398 * as well as the preferred amount for good performance.
399 * NB refers to the optimal block size for the immediately
400 * following subroutine, as returned by ILAENV.
401 *
402 MINWRK = 1
403 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
404 MINWRK = MAX( 10*( NMAX+1 ), 3*NMAX*NMAX )
405 NB = MAX( 1, ILAENV( 1, 'DGEQRF', ' ', NMAX, NMAX, -1, -1 ),
406 $ ILAENV( 1, 'DORMQR', 'LT', NMAX, NMAX, NMAX, -1 ),
407 $ ILAENV( 1, 'DORGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
408 MAXWRK = MAX( 10*( NMAX+1 ), 2*NMAX+NMAX*NB, 3*NMAX*NMAX )
409 WORK( 1 ) = MAXWRK
410 END IF
411 *
412 IF( LWORK.LT.MINWRK )
413 $ INFO = -20
414 *
415 IF( INFO.NE.0 ) THEN
416 CALL XERBLA( 'DDRGES', -INFO )
417 RETURN
418 END IF
419 *
420 * Quick return if possible
421 *
422 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
423 $ RETURN
424 *
425 SAFMIN = DLAMCH( 'Safe minimum' )
426 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
427 SAFMIN = SAFMIN / ULP
428 SAFMAX = ONE / SAFMIN
429 CALL DLABAD( SAFMIN, SAFMAX )
430 ULPINV = ONE / ULP
431 *
432 * The values RMAGN(2:3) depend on N, see below.
433 *
434 RMAGN( 0 ) = ZERO
435 RMAGN( 1 ) = ONE
436 *
437 * Loop over matrix sizes
438 *
439 NTESTT = 0
440 NERRS = 0
441 NMATS = 0
442 *
443 DO 190 JSIZE = 1, NSIZES
444 N = NN( JSIZE )
445 N1 = MAX( 1, N )
446 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
447 RMAGN( 3 ) = SAFMIN*ULPINV*DBLE( N1 )
448 *
449 IF( NSIZES.NE.1 ) THEN
450 MTYPES = MIN( MAXTYP, NTYPES )
451 ELSE
452 MTYPES = MIN( MAXTYP+1, NTYPES )
453 END IF
454 *
455 * Loop over matrix types
456 *
457 DO 180 JTYPE = 1, MTYPES
458 IF( .NOT.DOTYPE( JTYPE ) )
459 $ GO TO 180
460 NMATS = NMATS + 1
461 NTEST = 0
462 *
463 * Save ISEED in case of an error.
464 *
465 DO 20 J = 1, 4
466 IOLDSD( J ) = ISEED( J )
467 20 CONTINUE
468 *
469 * Initialize RESULT
470 *
471 DO 30 J = 1, 13
472 RESULT( J ) = ZERO
473 30 CONTINUE
474 *
475 * Generate test matrices A and B
476 *
477 * Description of control parameters:
478 *
479 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
480 * =3 means random.
481 * KATYPE: the "type" to be passed to DLATM4 for computing A.
482 * KAZERO: the pattern of zeros on the diagonal for A:
483 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
484 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
485 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
486 * non-zero entries.)
487 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
488 * =2: large, =3: small.
489 * IASIGN: 1 if the diagonal elements of A are to be
490 * multiplied by a random magnitude 1 number, =2 if
491 * randomly chosen diagonal blocks are to be rotated
492 * to form 2x2 blocks.
493 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
494 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
495 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
496 * RMAGN: used to implement KAMAGN and KBMAGN.
497 *
498 IF( MTYPES.GT.MAXTYP )
499 $ GO TO 110
500 IINFO = 0
501 IF( KCLASS( JTYPE ).LT.3 ) THEN
502 *
503 * Generate A (w/o rotation)
504 *
505 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
506 IN = 2*( ( N-1 ) / 2 ) + 1
507 IF( IN.NE.N )
508 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
509 ELSE
510 IN = N
511 END IF
512 CALL DLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
513 $ KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ),
514 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
515 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
516 $ ISEED, A, LDA )
517 IADD = KADD( KAZERO( JTYPE ) )
518 IF( IADD.GT.0 .AND. IADD.LE.N )
519 $ A( IADD, IADD ) = ONE
520 *
521 * Generate B (w/o rotation)
522 *
523 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
524 IN = 2*( ( N-1 ) / 2 ) + 1
525 IF( IN.NE.N )
526 $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDA )
527 ELSE
528 IN = N
529 END IF
530 CALL DLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
531 $ KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ),
532 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
533 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
534 $ ISEED, B, LDA )
535 IADD = KADD( KBZERO( JTYPE ) )
536 IF( IADD.NE.0 .AND. IADD.LE.N )
537 $ B( IADD, IADD ) = ONE
538 *
539 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
540 *
541 * Include rotations
542 *
543 * Generate Q, Z as Householder transformations times
544 * a diagonal matrix.
545 *
546 DO 50 JC = 1, N - 1
547 DO 40 JR = JC, N
548 Q( JR, JC ) = DLARND( 3, ISEED )
549 Z( JR, JC ) = DLARND( 3, ISEED )
550 40 CONTINUE
551 CALL DLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
552 $ WORK( JC ) )
553 WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) )
554 Q( JC, JC ) = ONE
555 CALL DLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
556 $ WORK( N+JC ) )
557 WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) )
558 Z( JC, JC ) = ONE
559 50 CONTINUE
560 Q( N, N ) = ONE
561 WORK( N ) = ZERO
562 WORK( 3*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
563 Z( N, N ) = ONE
564 WORK( 2*N ) = ZERO
565 WORK( 4*N ) = SIGN( ONE, DLARND( 2, ISEED ) )
566 *
567 * Apply the diagonal matrices
568 *
569 DO 70 JC = 1, N
570 DO 60 JR = 1, N
571 A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
572 $ A( JR, JC )
573 B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
574 $ B( JR, JC )
575 60 CONTINUE
576 70 CONTINUE
577 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
578 $ LDA, WORK( 2*N+1 ), IINFO )
579 IF( IINFO.NE.0 )
580 $ GO TO 100
581 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
582 $ A, LDA, WORK( 2*N+1 ), IINFO )
583 IF( IINFO.NE.0 )
584 $ GO TO 100
585 CALL DORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
586 $ LDA, WORK( 2*N+1 ), IINFO )
587 IF( IINFO.NE.0 )
588 $ GO TO 100
589 CALL DORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
590 $ B, LDA, WORK( 2*N+1 ), IINFO )
591 IF( IINFO.NE.0 )
592 $ GO TO 100
593 END IF
594 ELSE
595 *
596 * Random matrices
597 *
598 DO 90 JC = 1, N
599 DO 80 JR = 1, N
600 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
601 $ DLARND( 2, ISEED )
602 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
603 $ DLARND( 2, ISEED )
604 80 CONTINUE
605 90 CONTINUE
606 END IF
607 *
608 100 CONTINUE
609 *
610 IF( IINFO.NE.0 ) THEN
611 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
612 $ IOLDSD
613 INFO = ABS( IINFO )
614 RETURN
615 END IF
616 *
617 110 CONTINUE
618 *
619 DO 120 I = 1, 13
620 RESULT( I ) = -ONE
621 120 CONTINUE
622 *
623 * Test with and without sorting of eigenvalues
624 *
625 DO 150 ISORT = 0, 1
626 IF( ISORT.EQ.0 ) THEN
627 SORT = 'N'
628 RSUB = 0
629 ELSE
630 SORT = 'S'
631 RSUB = 5
632 END IF
633 *
634 * Call DGGES to compute H, T, Q, Z, alpha, and beta.
635 *
636 CALL DLACPY( 'Full', N, N, A, LDA, S, LDA )
637 CALL DLACPY( 'Full', N, N, B, LDA, T, LDA )
638 NTEST = 1 + RSUB + ISORT
639 RESULT( 1+RSUB+ISORT ) = ULPINV
640 CALL DGGES( 'V', 'V', SORT, DLCTES, N, S, LDA, T, LDA,
641 $ SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDQ,
642 $ WORK, LWORK, BWORK, IINFO )
643 IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
644 RESULT( 1+RSUB+ISORT ) = ULPINV
645 WRITE( NOUNIT, FMT = 9999 )'DGGES', IINFO, N, JTYPE,
646 $ IOLDSD
647 INFO = ABS( IINFO )
648 GO TO 160
649 END IF
650 *
651 NTEST = 4 + RSUB
652 *
653 * Do tests 1--4 (or tests 7--9 when reordering )
654 *
655 IF( ISORT.EQ.0 ) THEN
656 CALL DGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ,
657 $ WORK, RESULT( 1 ) )
658 CALL DGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ,
659 $ WORK, RESULT( 2 ) )
660 ELSE
661 CALL DGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q,
662 $ LDQ, Z, LDQ, WORK, RESULT( 7 ) )
663 END IF
664 CALL DGET51( 3, N, A, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
665 $ RESULT( 3+RSUB ) )
666 CALL DGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
667 $ RESULT( 4+RSUB ) )
668 *
669 * Do test 5 and 6 (or Tests 10 and 11 when reordering):
670 * check Schur form of A and compare eigenvalues with
671 * diagonals.
672 *
673 NTEST = 6 + RSUB
674 TEMP1 = ZERO
675 *
676 DO 130 J = 1, N
677 ILABAD = .FALSE.
678 IF( ALPHAI( J ).EQ.ZERO ) THEN
679 TEMP2 = ( ABS( ALPHAR( J )-S( J, J ) ) /
680 $ MAX( SAFMIN, ABS( ALPHAR( J ) ), ABS( S( J,
681 $ J ) ) )+ABS( BETA( J )-T( J, J ) ) /
682 $ MAX( SAFMIN, ABS( BETA( J ) ), ABS( T( J,
683 $ J ) ) ) ) / ULP
684 *
685 IF( J.LT.N ) THEN
686 IF( S( J+1, J ).NE.ZERO ) THEN
687 ILABAD = .TRUE.
688 RESULT( 5+RSUB ) = ULPINV
689 END IF
690 END IF
691 IF( J.GT.1 ) THEN
692 IF( S( J, J-1 ).NE.ZERO ) THEN
693 ILABAD = .TRUE.
694 RESULT( 5+RSUB ) = ULPINV
695 END IF
696 END IF
697 *
698 ELSE
699 IF( ALPHAI( J ).GT.ZERO ) THEN
700 I1 = J
701 ELSE
702 I1 = J - 1
703 END IF
704 IF( I1.LE.0 .OR. I1.GE.N ) THEN
705 ILABAD = .TRUE.
706 ELSE IF( I1.LT.N-1 ) THEN
707 IF( S( I1+2, I1+1 ).NE.ZERO ) THEN
708 ILABAD = .TRUE.
709 RESULT( 5+RSUB ) = ULPINV
710 END IF
711 ELSE IF( I1.GT.1 ) THEN
712 IF( S( I1, I1-1 ).NE.ZERO ) THEN
713 ILABAD = .TRUE.
714 RESULT( 5+RSUB ) = ULPINV
715 END IF
716 END IF
717 IF( .NOT.ILABAD ) THEN
718 CALL DGET53( S( I1, I1 ), LDA, T( I1, I1 ), LDA,
719 $ BETA( J ), ALPHAR( J ),
720 $ ALPHAI( J ), TEMP2, IERR )
721 IF( IERR.GE.3 ) THEN
722 WRITE( NOUNIT, FMT = 9998 )IERR, J, N,
723 $ JTYPE, IOLDSD
724 INFO = ABS( IERR )
725 END IF
726 ELSE
727 TEMP2 = ULPINV
728 END IF
729 *
730 END IF
731 TEMP1 = MAX( TEMP1, TEMP2 )
732 IF( ILABAD ) THEN
733 WRITE( NOUNIT, FMT = 9997 )J, N, JTYPE, IOLDSD
734 END IF
735 130 CONTINUE
736 RESULT( 6+RSUB ) = TEMP1
737 *
738 IF( ISORT.GE.1 ) THEN
739 *
740 * Do test 12
741 *
742 NTEST = 12
743 RESULT( 12 ) = ZERO
744 KNTEIG = 0
745 DO 140 I = 1, N
746 IF( DLCTES( ALPHAR( I ), ALPHAI( I ),
747 $ BETA( I ) ) .OR. DLCTES( ALPHAR( I ),
748 $ -ALPHAI( I ), BETA( I ) ) ) THEN
749 KNTEIG = KNTEIG + 1
750 END IF
751 IF( I.LT.N ) THEN
752 IF( ( DLCTES( ALPHAR( I+1 ), ALPHAI( I+1 ),
753 $ BETA( I+1 ) ) .OR. DLCTES( ALPHAR( I+1 ),
754 $ -ALPHAI( I+1 ), BETA( I+1 ) ) ) .AND.
755 $ ( .NOT.( DLCTES( ALPHAR( I ), ALPHAI( I ),
756 $ BETA( I ) ) .OR. DLCTES( ALPHAR( I ),
757 $ -ALPHAI( I ), BETA( I ) ) ) ) .AND.
758 $ IINFO.NE.N+2 ) THEN
759 RESULT( 12 ) = ULPINV
760 END IF
761 END IF
762 140 CONTINUE
763 IF( SDIM.NE.KNTEIG ) THEN
764 RESULT( 12 ) = ULPINV
765 END IF
766 END IF
767 *
768 150 CONTINUE
769 *
770 * End of Loop -- Check for RESULT(j) > THRESH
771 *
772 160 CONTINUE
773 *
774 NTESTT = NTESTT + NTEST
775 *
776 * Print out tests which fail.
777 *
778 DO 170 JR = 1, NTEST
779 IF( RESULT( JR ).GE.THRESH ) THEN
780 *
781 * If this is the first test to fail,
782 * print a header to the data file.
783 *
784 IF( NERRS.EQ.0 ) THEN
785 WRITE( NOUNIT, FMT = 9996 )'DGS'
786 *
787 * Matrix types
788 *
789 WRITE( NOUNIT, FMT = 9995 )
790 WRITE( NOUNIT, FMT = 9994 )
791 WRITE( NOUNIT, FMT = 9993 )'Orthogonal'
792 *
793 * Tests performed
794 *
795 WRITE( NOUNIT, FMT = 9992 )'orthogonal', '''',
796 $ 'transpose', ( '''', J = 1, 8 )
797 *
798 END IF
799 NERRS = NERRS + 1
800 IF( RESULT( JR ).LT.10000.0D0 ) THEN
801 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
802 $ RESULT( JR )
803 ELSE
804 WRITE( NOUNIT, FMT = 9990 )N, JTYPE, IOLDSD, JR,
805 $ RESULT( JR )
806 END IF
807 END IF
808 170 CONTINUE
809 *
810 180 CONTINUE
811 190 CONTINUE
812 *
813 * Summary
814 *
815 CALL ALASVM( 'DGS', NOUNIT, NERRS, NTESTT, 0 )
816 *
817 WORK( 1 ) = MAXWRK
818 *
819 RETURN
820 *
821 9999 FORMAT( ' DDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
822 $ I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
823 *
824 9998 FORMAT( ' DDRGES: DGET53 returned INFO=', I1, ' for eigenvalue ',
825 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(',
826 $ 4( I4, ',' ), I5, ')' )
827 *
828 9997 FORMAT( ' DDRGES: S not in Schur form at eigenvalue ', I6, '.',
829 $ / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
830 $ I5, ')' )
831 *
832 9996 FORMAT( / 1X, A3, ' -- Real Generalized Schur form driver' )
833 *
834 9995 FORMAT( ' Matrix types (see DDRGES for details): ' )
835 *
836 9994 FORMAT( ' Special Matrices:', 23X,
837 $ '(J''=transposed Jordan block)',
838 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
839 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
840 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
841 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
842 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
843 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
844 9993 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
845 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
846 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
847 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
848 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
849 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
850 $ '23=(small,large) 24=(small,small) 25=(large,large)',
851 $ / ' 26=random O(1) matrices.' )
852 *
853 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
854 $ 'Q and Z are ', A, ',', / 19X,
855 $ 'l and r are the appropriate left and right', / 19X,
856 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
857 $ ' means ', A, '.)', / ' Without ordering: ',
858 $ / ' 1 = | A - Q S Z', A,
859 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
860 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
861 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
862 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
863 $ / ' 6 = difference between (alpha,beta)',
864 $ ' and diagonals of (S,T)', / ' With ordering: ',
865 $ / ' 7 = | (A,B) - Q (S,T) Z', A,
866 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', A,
867 $ ' | / ( n ulp ) 9 = | I - ZZ', A,
868 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
869 $ / ' 11 = difference between (alpha,beta) and diagonals',
870 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
871 $ 'selected eigenvalues', / )
872 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
873 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
874 9990 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
875 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
876 *
877 * End of DDRGES
878 *
879 END