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