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