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