1 SUBROUTINE SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
2 $ AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
3 $ WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
11 $ NOUT, NSIZE
12 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 INTEGER IWORK( * )
17 REAL A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
18 $ ALPHAR( * ), B( LDA, * ), BETA( * ),
19 $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
20 $ WORK( * ), Z( LDA, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem expert driver SGGESX.
28 *
29 * SGGESX factors A and B as Q S Z' and Q T Z', where ' means
30 * transpose, T is upper triangular, S is in generalized Schur form
31 * (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
32 * the 2x2 blocks corresponding to complex conjugate pairs of
33 * generalized eigenvalues), and Q and Z are orthogonal. It also
34 * computes the generalized eigenvalues (alpha(1),beta(1)), ...,
35 * (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
36 * characteristic equation
37 *
38 * det( A - w(j) B ) = 0
39 *
40 * Optionally it also reorders the eigenvalues so that a selected
41 * cluster of eigenvalues appears in the leading diagonal block of the
42 * Schur forms; computes a reciprocal condition number for the average
43 * of the selected eigenvalues; and computes a reciprocal condition
44 * number for the right and left deflating subspaces corresponding to
45 * the selected eigenvalues.
46 *
47 * When SDRGSX is called with NSIZE > 0, five (5) types of built-in
48 * matrix pairs are used to test the routine SGGESX.
49 *
50 * When SDRGSX is called with NSIZE = 0, it reads in test matrix data
51 * to test SGGESX.
52 *
53 * For each matrix pair, the following tests will be performed and
54 * compared with the threshhold THRESH except for the tests (7) and (9):
55 *
56 * (1) | A - Q S Z' | / ( |A| n ulp )
57 *
58 * (2) | B - Q T Z' | / ( |B| n ulp )
59 *
60 * (3) | I - QQ' | / ( n ulp )
61 *
62 * (4) | I - ZZ' | / ( n ulp )
63 *
64 * (5) if A is in Schur form (i.e. quasi-triangular form)
65 *
66 * (6) maximum over j of D(j) where:
67 *
68 * if alpha(j) is real:
69 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
70 * D(j) = ------------------------ + -----------------------
71 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
72 *
73 * if alpha(j) is complex:
74 * | det( s S - w T ) |
75 * D(j) = ---------------------------------------------------
76 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
77 *
78 * and S and T are here the 2 x 2 diagonal blocks of S and T
79 * corresponding to the j-th and j+1-th eigenvalues.
80 *
81 * (7) if sorting worked and SDIM is the number of eigenvalues
82 * which were selected.
83 *
84 * (8) the estimated value DIF does not differ from the true values of
85 * Difu and Difl more than a factor 10*THRESH. If the estimate DIF
86 * equals zero the corresponding true values of Difu and Difl
87 * should be less than EPS*norm(A, B). If the true value of Difu
88 * and Difl equal zero, the estimate DIF should be less than
89 * EPS*norm(A, B).
90 *
91 * (9) If INFO = N+3 is returned by SGGESX, the reordering "failed"
92 * and we check that DIF = PL = PR = 0 and that the true value of
93 * Difu and Difl is < EPS*norm(A, B). We count the events when
94 * INFO=N+3.
95 *
96 * For read-in test matrices, the above tests are run except that the
97 * exact value for DIF (and PL) is input data. Additionally, there is
98 * one more test run for read-in test matrices:
99 *
100 * (10) the estimated value PL does not differ from the true value of
101 * PLTRU more than a factor THRESH. If the estimate PL equals
102 * zero the corresponding true value of PLTRU should be less than
103 * EPS*norm(A, B). If the true value of PLTRU equal zero, the
104 * estimate PL should be less than EPS*norm(A, B).
105 *
106 * Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
107 * matrix pairs are generated and tested. NSIZE should be kept small.
108 *
109 * SVD (routine SGESVD) is used for computing the true value of DIF_u
110 * and DIF_l when testing the built-in test problems.
111 *
112 * Built-in Test Matrices
113 * ======================
114 *
115 * All built-in test matrices are the 2 by 2 block of triangular
116 * matrices
117 *
118 * A = [ A11 A12 ] and B = [ B11 B12 ]
119 * [ A22 ] [ B22 ]
120 *
121 * where for different type of A11 and A22 are given as the following.
122 * A12 and B12 are chosen so that the generalized Sylvester equation
123 *
124 * A11*R - L*A22 = -A12
125 * B11*R - L*B22 = -B12
126 *
127 * have prescribed solution R and L.
128 *
129 * Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
130 * B11 = I_m, B22 = I_k
131 * where J_k(a,b) is the k-by-k Jordan block with ``a'' on
132 * diagonal and ``b'' on superdiagonal.
133 *
134 * Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and
135 * B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
136 * A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
137 * B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
138 *
139 * Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each
140 * second diagonal block in A_11 and each third diagonal block
141 * in A_22 are made as 2 by 2 blocks.
142 *
143 * Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
144 * for i=1,...,m, j=1,...,m and
145 * A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
146 * for i=m+1,...,k, j=m+1,...,k
147 *
148 * Type 5: (A,B) and have potentially close or common eigenvalues and
149 * very large departure from block diagonality A_11 is chosen
150 * as the m x m leading submatrix of A_1:
151 * | 1 b |
152 * | -b 1 |
153 * | 1+d b |
154 * | -b 1+d |
155 * A_1 = | d 1 |
156 * | -1 d |
157 * | -d 1 |
158 * | -1 -d |
159 * | 1 |
160 * and A_22 is chosen as the k x k leading submatrix of A_2:
161 * | -1 b |
162 * | -b -1 |
163 * | 1-d b |
164 * | -b 1-d |
165 * A_2 = | d 1+b |
166 * | -1-b d |
167 * | -d 1+b |
168 * | -1+b -d |
169 * | 1-d |
170 * and matrix B are chosen as identity matrices (see SLATM5).
171 *
172 *
173 * Arguments
174 * =========
175 *
176 * NSIZE (input) INTEGER
177 * The maximum size of the matrices to use. NSIZE >= 0.
178 * If NSIZE = 0, no built-in tests matrices are used, but
179 * read-in test matrices are used to test SGGESX.
180 *
181 * NCMAX (input) INTEGER
182 * Maximum allowable NMAX for generating Kroneker matrix
183 * in call to SLAKF2
184 *
185 * THRESH (input) REAL
186 * A test will count as "failed" if the "error", computed as
187 * described above, exceeds THRESH. Note that the error
188 * is scaled to be O(1), so THRESH should be a reasonably
189 * small multiple of 1, e.g., 10 or 100. In particular,
190 * it should not depend on the precision (single vs. double)
191 * or the size of the matrix. THRESH >= 0.
192 *
193 * NIN (input) INTEGER
194 * The FORTRAN unit number for reading in the data file of
195 * problems to solve.
196 *
197 * NOUT (input) INTEGER
198 * The FORTRAN unit number for printing out error messages
199 * (e.g., if a routine returns IINFO not equal to 0.)
200 *
201 * A (workspace) REAL array, dimension (LDA, NSIZE)
202 * Used to store the matrix whose eigenvalues are to be
203 * computed. On exit, A contains the last matrix actually used.
204 *
205 * LDA (input) INTEGER
206 * The leading dimension of A, B, AI, BI, Z and Q,
207 * LDA >= max( 1, NSIZE ). For the read-in test,
208 * LDA >= max( 1, N ), N is the size of the test matrices.
209 *
210 * B (workspace) REAL array, dimension (LDA, NSIZE)
211 * Used to store the matrix whose eigenvalues are to be
212 * computed. On exit, B contains the last matrix actually used.
213 *
214 * AI (workspace) REAL array, dimension (LDA, NSIZE)
215 * Copy of A, modified by SGGESX.
216 *
217 * BI (workspace) REAL array, dimension (LDA, NSIZE)
218 * Copy of B, modified by SGGESX.
219 *
220 * Z (workspace) REAL array, dimension (LDA, NSIZE)
221 * Z holds the left Schur vectors computed by SGGESX.
222 *
223 * Q (workspace) REAL array, dimension (LDA, NSIZE)
224 * Q holds the right Schur vectors computed by SGGESX.
225 *
226 * ALPHAR (workspace) REAL array, dimension (NSIZE)
227 * ALPHAI (workspace) REAL array, dimension (NSIZE)
228 * BETA (workspace) REAL array, dimension (NSIZE)
229 * On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
230 *
231 * C (workspace) REAL array, dimension (LDC, LDC)
232 * Store the matrix generated by subroutine SLAKF2, this is the
233 * matrix formed by Kronecker products used for estimating
234 * DIF.
235 *
236 * LDC (input) INTEGER
237 * The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
238 *
239 * S (workspace) REAL array, dimension (LDC)
240 * Singular values of C
241 *
242 * WORK (workspace) REAL array, dimension (LWORK)
243 *
244 * LWORK (input) INTEGER
245 * The dimension of the array WORK.
246 * LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
247 *
248 * IWORK (workspace) INTEGER array, dimension (LIWORK)
249 *
250 * LIWORK (input) INTEGER
251 * The dimension of the array IWORK. LIWORK >= NSIZE + 6.
252 *
253 * BWORK (workspace) LOGICAL array, dimension (LDA)
254 *
255 * INFO (output) INTEGER
256 * = 0: successful exit
257 * < 0: if INFO = -i, the i-th argument had an illegal value.
258 * > 0: A routine returned an error code.
259 *
260 * =====================================================================
261 *
262 * .. Parameters ..
263 REAL ZERO, ONE, TEN
264 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
265 * ..
266 * .. Local Scalars ..
267 LOGICAL ILABAD
268 CHARACTER SENSE
269 INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
270 $ MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
271 $ PRTYPE, QBA, QBB
272 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
273 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
274 * ..
275 * .. Local Arrays ..
276 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
277 * ..
278 * .. External Functions ..
279 LOGICAL SLCTSX
280 INTEGER ILAENV
281 REAL SLAMCH, SLANGE
282 EXTERNAL SLCTSX, ILAENV, SLAMCH, SLANGE
283 * ..
284 * .. External Subroutines ..
285 EXTERNAL ALASVM, SGESVD, SGET51, SGET53, SGGESX, SLABAD,
286 $ SLACPY, SLAKF2, SLASET, SLATM5, XERBLA
287 * ..
288 * .. Intrinsic Functions ..
289 INTRINSIC ABS, MAX, SQRT
290 * ..
291 * .. Scalars in Common ..
292 LOGICAL FS
293 INTEGER K, M, MPLUSN, N
294 * ..
295 * .. Common blocks ..
296 COMMON / MN / M, N, MPLUSN, K, FS
297 * ..
298 * .. Executable Statements ..
299 *
300 * Check for errors
301 *
302 IF( NSIZE.LT.0 ) THEN
303 INFO = -1
304 ELSE IF( THRESH.LT.ZERO ) THEN
305 INFO = -2
306 ELSE IF( NIN.LE.0 ) THEN
307 INFO = -3
308 ELSE IF( NOUT.LE.0 ) THEN
309 INFO = -4
310 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
311 INFO = -6
312 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
313 INFO = -17
314 ELSE IF( LIWORK.LT.NSIZE+6 ) THEN
315 INFO = -21
316 END IF
317 *
318 * Compute workspace
319 * (Note: Comments in the code beginning "Workspace:" describe the
320 * minimal amount of workspace needed at that point in the code,
321 * as well as the preferred amount for good performance.
322 * NB refers to the optimal block size for the immediately
323 * following subroutine, as returned by ILAENV.)
324 *
325 MINWRK = 1
326 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
327 c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
328 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
329 *
330 * workspace for sggesx
331 *
332 MAXWRK = 9*( NSIZE+1 ) + NSIZE*
333 $ ILAENV( 1, 'SGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
334 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
335 $ ILAENV( 1, 'SORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
336 *
337 * workspace for sgesvd
338 *
339 BDSPAC = 5*NSIZE*NSIZE / 2
340 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
341 $ ILAENV( 1, 'SGEBRD', ' ', NSIZE*NSIZE / 2,
342 $ NSIZE*NSIZE / 2, -1, -1 ) )
343 MAXWRK = MAX( MAXWRK, BDSPAC )
344 *
345 MAXWRK = MAX( MAXWRK, MINWRK )
346 *
347 WORK( 1 ) = MAXWRK
348 END IF
349 *
350 IF( LWORK.LT.MINWRK )
351 $ INFO = -19
352 *
353 IF( INFO.NE.0 ) THEN
354 CALL XERBLA( 'SDRGSX', -INFO )
355 RETURN
356 END IF
357 *
358 * Important constants
359 *
360 ULP = SLAMCH( 'P' )
361 ULPINV = ONE / ULP
362 SMLNUM = SLAMCH( 'S' ) / ULP
363 BIGNUM = ONE / SMLNUM
364 CALL SLABAD( SMLNUM, BIGNUM )
365 THRSH2 = TEN*THRESH
366 NTESTT = 0
367 NERRS = 0
368 *
369 * Go to the tests for read-in matrix pairs
370 *
371 IFUNC = 0
372 IF( NSIZE.EQ.0 )
373 $ GO TO 70
374 *
375 * Test the built-in matrix pairs.
376 * Loop over different functions (IFUNC) of SGGESX, types (PRTYPE)
377 * of test matrices, different size (M+N)
378 *
379 PRTYPE = 0
380 QBA = 3
381 QBB = 4
382 WEIGHT = SQRT( ULP )
383 *
384 DO 60 IFUNC = 0, 3
385 DO 50 PRTYPE = 1, 5
386 DO 40 M = 1, NSIZE - 1
387 DO 30 N = 1, NSIZE - M
388 *
389 WEIGHT = ONE / WEIGHT
390 MPLUSN = M + N
391 *
392 * Generate test matrices
393 *
394 FS = .TRUE.
395 K = 0
396 *
397 CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
398 $ LDA )
399 CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
400 $ LDA )
401 *
402 CALL SLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
403 $ LDA, AI( 1, M+1 ), LDA, BI, LDA,
404 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
405 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
406 *
407 * Compute the Schur factorization and swapping the
408 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
409 * Swapping is accomplished via the function SLCTSX
410 * which is supplied below.
411 *
412 IF( IFUNC.EQ.0 ) THEN
413 SENSE = 'N'
414 ELSE IF( IFUNC.EQ.1 ) THEN
415 SENSE = 'E'
416 ELSE IF( IFUNC.EQ.2 ) THEN
417 SENSE = 'V'
418 ELSE IF( IFUNC.EQ.3 ) THEN
419 SENSE = 'B'
420 END IF
421 *
422 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
423 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
424 *
425 CALL SGGESX( 'V', 'V', 'S', SLCTSX, SENSE, MPLUSN, AI,
426 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
427 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
428 $ IWORK, LIWORK, BWORK, LINFO )
429 *
430 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
431 RESULT( 1 ) = ULPINV
432 WRITE( NOUT, FMT = 9999 )'SGGESX', LINFO, MPLUSN,
433 $ PRTYPE
434 INFO = LINFO
435 GO TO 30
436 END IF
437 *
438 * Compute the norm(A, B)
439 *
440 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
441 $ MPLUSN )
442 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
443 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
444 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
445 $ WORK )
446 *
447 * Do tests (1) to (4)
448 *
449 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
450 $ LDA, WORK, RESULT( 1 ) )
451 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
452 $ LDA, WORK, RESULT( 2 ) )
453 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
454 $ LDA, WORK, RESULT( 3 ) )
455 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
456 $ LDA, WORK, RESULT( 4 ) )
457 NTEST = 4
458 *
459 * Do tests (5) and (6): check Schur form of A and
460 * compare eigenvalues with diagonals.
461 *
462 TEMP1 = ZERO
463 RESULT( 5 ) = ZERO
464 RESULT( 6 ) = ZERO
465 *
466 DO 10 J = 1, MPLUSN
467 ILABAD = .FALSE.
468 IF( ALPHAI( J ).EQ.ZERO ) THEN
469 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
470 $ MAX( SMLNUM, ABS( ALPHAR( J ) ),
471 $ ABS( AI( J, J ) ) )+
472 $ ABS( BETA( J )-BI( J, J ) ) /
473 $ MAX( SMLNUM, ABS( BETA( J ) ),
474 $ ABS( BI( J, J ) ) ) ) / ULP
475 IF( J.LT.MPLUSN ) THEN
476 IF( AI( J+1, J ).NE.ZERO ) THEN
477 ILABAD = .TRUE.
478 RESULT( 5 ) = ULPINV
479 END IF
480 END IF
481 IF( J.GT.1 ) THEN
482 IF( AI( J, J-1 ).NE.ZERO ) THEN
483 ILABAD = .TRUE.
484 RESULT( 5 ) = ULPINV
485 END IF
486 END IF
487 ELSE
488 IF( ALPHAI( J ).GT.ZERO ) THEN
489 I1 = J
490 ELSE
491 I1 = J - 1
492 END IF
493 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
494 ILABAD = .TRUE.
495 ELSE IF( I1.LT.MPLUSN-1 ) THEN
496 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
497 ILABAD = .TRUE.
498 RESULT( 5 ) = ULPINV
499 END IF
500 ELSE IF( I1.GT.1 ) THEN
501 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
502 ILABAD = .TRUE.
503 RESULT( 5 ) = ULPINV
504 END IF
505 END IF
506 IF( .NOT.ILABAD ) THEN
507 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
508 $ LDA, BETA( J ), ALPHAR( J ),
509 $ ALPHAI( J ), TEMP2, IINFO )
510 IF( IINFO.GE.3 ) THEN
511 WRITE( NOUT, FMT = 9997 )IINFO, J,
512 $ MPLUSN, PRTYPE
513 INFO = ABS( IINFO )
514 END IF
515 ELSE
516 TEMP2 = ULPINV
517 END IF
518 END IF
519 TEMP1 = MAX( TEMP1, TEMP2 )
520 IF( ILABAD ) THEN
521 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
522 END IF
523 10 CONTINUE
524 RESULT( 6 ) = TEMP1
525 NTEST = NTEST + 2
526 *
527 * Test (7) (if sorting worked)
528 *
529 RESULT( 7 ) = ZERO
530 IF( LINFO.EQ.MPLUSN+3 ) THEN
531 RESULT( 7 ) = ULPINV
532 ELSE IF( MM.NE.N ) THEN
533 RESULT( 7 ) = ULPINV
534 END IF
535 NTEST = NTEST + 1
536 *
537 * Test (8): compare the estimated value DIF and its
538 * value. first, compute the exact DIF.
539 *
540 RESULT( 8 ) = ZERO
541 MN2 = MM*( MPLUSN-MM )*2
542 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
543 *
544 * Note: for either following two causes, there are
545 * almost same number of test cases fail the test.
546 *
547 CALL SLAKF2( MM, MPLUSN-MM, AI, LDA,
548 $ AI( MM+1, MM+1 ), BI,
549 $ BI( MM+1, MM+1 ), C, LDC )
550 *
551 CALL SGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
552 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
553 $ INFO )
554 DIFTRU = S( MN2 )
555 *
556 IF( DIFEST( 2 ).EQ.ZERO ) THEN
557 IF( DIFTRU.GT.ABNRM*ULP )
558 $ RESULT( 8 ) = ULPINV
559 ELSE IF( DIFTRU.EQ.ZERO ) THEN
560 IF( DIFEST( 2 ).GT.ABNRM*ULP )
561 $ RESULT( 8 ) = ULPINV
562 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
563 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
564 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
565 $ DIFEST( 2 ) / DIFTRU )
566 END IF
567 NTEST = NTEST + 1
568 END IF
569 *
570 * Test (9)
571 *
572 RESULT( 9 ) = ZERO
573 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
574 IF( DIFTRU.GT.ABNRM*ULP )
575 $ RESULT( 9 ) = ULPINV
576 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
577 $ RESULT( 9 ) = ULPINV
578 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
579 $ RESULT( 9 ) = ULPINV
580 NTEST = NTEST + 1
581 END IF
582 *
583 NTESTT = NTESTT + NTEST
584 *
585 * Print out tests which fail.
586 *
587 DO 20 J = 1, 9
588 IF( RESULT( J ).GE.THRESH ) THEN
589 *
590 * If this is the first test to fail,
591 * print a header to the data file.
592 *
593 IF( NERRS.EQ.0 ) THEN
594 WRITE( NOUT, FMT = 9995 )'SGX'
595 *
596 * Matrix types
597 *
598 WRITE( NOUT, FMT = 9993 )
599 *
600 * Tests performed
601 *
602 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
603 $ 'transpose', ( '''', I = 1, 4 )
604 *
605 END IF
606 NERRS = NERRS + 1
607 IF( RESULT( J ).LT.10000.0 ) THEN
608 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
609 $ WEIGHT, M, J, RESULT( J )
610 ELSE
611 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
612 $ WEIGHT, M, J, RESULT( J )
613 END IF
614 END IF
615 20 CONTINUE
616 *
617 30 CONTINUE
618 40 CONTINUE
619 50 CONTINUE
620 60 CONTINUE
621 *
622 GO TO 150
623 *
624 70 CONTINUE
625 *
626 * Read in data from file to check accuracy of condition estimation
627 * Read input data until N=0
628 *
629 NPTKNT = 0
630 *
631 80 CONTINUE
632 READ( NIN, FMT = *, END = 140 )MPLUSN
633 IF( MPLUSN.EQ.0 )
634 $ GO TO 140
635 READ( NIN, FMT = *, END = 140 )N
636 DO 90 I = 1, MPLUSN
637 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
638 90 CONTINUE
639 DO 100 I = 1, MPLUSN
640 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
641 100 CONTINUE
642 READ( NIN, FMT = * )PLTRU, DIFTRU
643 *
644 NPTKNT = NPTKNT + 1
645 FS = .TRUE.
646 K = 0
647 M = MPLUSN - N
648 *
649 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
650 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
651 *
652 * Compute the Schur factorization while swaping the
653 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
654 *
655 CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
656 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
657 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
658 *
659 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
660 RESULT( 1 ) = ULPINV
661 WRITE( NOUT, FMT = 9998 )'SGGESX', LINFO, MPLUSN, NPTKNT
662 GO TO 130
663 END IF
664 *
665 * Compute the norm(A, B)
666 * (should this be norm of (A,B) or (AI,BI)?)
667 *
668 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
669 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
670 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
671 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
672 *
673 * Do tests (1) to (4)
674 *
675 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
676 $ RESULT( 1 ) )
677 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
678 $ RESULT( 2 ) )
679 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
680 $ RESULT( 3 ) )
681 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
682 $ RESULT( 4 ) )
683 *
684 * Do tests (5) and (6): check Schur form of A and compare
685 * eigenvalues with diagonals.
686 *
687 NTEST = 6
688 TEMP1 = ZERO
689 RESULT( 5 ) = ZERO
690 RESULT( 6 ) = ZERO
691 *
692 DO 110 J = 1, MPLUSN
693 ILABAD = .FALSE.
694 IF( ALPHAI( J ).EQ.ZERO ) THEN
695 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
696 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
697 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
698 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
699 $ / ULP
700 IF( J.LT.MPLUSN ) THEN
701 IF( AI( J+1, J ).NE.ZERO ) THEN
702 ILABAD = .TRUE.
703 RESULT( 5 ) = ULPINV
704 END IF
705 END IF
706 IF( J.GT.1 ) THEN
707 IF( AI( J, J-1 ).NE.ZERO ) THEN
708 ILABAD = .TRUE.
709 RESULT( 5 ) = ULPINV
710 END IF
711 END IF
712 ELSE
713 IF( ALPHAI( J ).GT.ZERO ) THEN
714 I1 = J
715 ELSE
716 I1 = J - 1
717 END IF
718 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
719 ILABAD = .TRUE.
720 ELSE IF( I1.LT.MPLUSN-1 ) THEN
721 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
722 ILABAD = .TRUE.
723 RESULT( 5 ) = ULPINV
724 END IF
725 ELSE IF( I1.GT.1 ) THEN
726 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
727 ILABAD = .TRUE.
728 RESULT( 5 ) = ULPINV
729 END IF
730 END IF
731 IF( .NOT.ILABAD ) THEN
732 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
733 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
734 $ IINFO )
735 IF( IINFO.GE.3 ) THEN
736 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
737 INFO = ABS( IINFO )
738 END IF
739 ELSE
740 TEMP2 = ULPINV
741 END IF
742 END IF
743 TEMP1 = MAX( TEMP1, TEMP2 )
744 IF( ILABAD ) THEN
745 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
746 END IF
747 110 CONTINUE
748 RESULT( 6 ) = TEMP1
749 *
750 * Test (7) (if sorting worked) <--------- need to be checked.
751 *
752 NTEST = 7
753 RESULT( 7 ) = ZERO
754 IF( LINFO.EQ.MPLUSN+3 )
755 $ RESULT( 7 ) = ULPINV
756 *
757 * Test (8): compare the estimated value of DIF and its true value.
758 *
759 NTEST = 8
760 RESULT( 8 ) = ZERO
761 IF( DIFEST( 2 ).EQ.ZERO ) THEN
762 IF( DIFTRU.GT.ABNRM*ULP )
763 $ RESULT( 8 ) = ULPINV
764 ELSE IF( DIFTRU.EQ.ZERO ) THEN
765 IF( DIFEST( 2 ).GT.ABNRM*ULP )
766 $ RESULT( 8 ) = ULPINV
767 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
768 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
769 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
770 END IF
771 *
772 * Test (9)
773 *
774 NTEST = 9
775 RESULT( 9 ) = ZERO
776 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
777 IF( DIFTRU.GT.ABNRM*ULP )
778 $ RESULT( 9 ) = ULPINV
779 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
780 $ RESULT( 9 ) = ULPINV
781 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
782 $ RESULT( 9 ) = ULPINV
783 END IF
784 *
785 * Test (10): compare the estimated value of PL and it true value.
786 *
787 NTEST = 10
788 RESULT( 10 ) = ZERO
789 IF( PL( 1 ).EQ.ZERO ) THEN
790 IF( PLTRU.GT.ABNRM*ULP )
791 $ RESULT( 10 ) = ULPINV
792 ELSE IF( PLTRU.EQ.ZERO ) THEN
793 IF( PL( 1 ).GT.ABNRM*ULP )
794 $ RESULT( 10 ) = ULPINV
795 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
796 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
797 RESULT( 10 ) = ULPINV
798 END IF
799 *
800 NTESTT = NTESTT + NTEST
801 *
802 * Print out tests which fail.
803 *
804 DO 120 J = 1, NTEST
805 IF( RESULT( J ).GE.THRESH ) THEN
806 *
807 * If this is the first test to fail,
808 * print a header to the data file.
809 *
810 IF( NERRS.EQ.0 ) THEN
811 WRITE( NOUT, FMT = 9995 )'SGX'
812 *
813 * Matrix types
814 *
815 WRITE( NOUT, FMT = 9994 )
816 *
817 * Tests performed
818 *
819 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
820 $ 'transpose', ( '''', I = 1, 4 )
821 *
822 END IF
823 NERRS = NERRS + 1
824 IF( RESULT( J ).LT.10000.0 ) THEN
825 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
826 ELSE
827 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
828 END IF
829 END IF
830 *
831 120 CONTINUE
832 *
833 130 CONTINUE
834 GO TO 80
835 140 CONTINUE
836 *
837 150 CONTINUE
838 *
839 * Summary
840 *
841 CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
842 *
843 WORK( 1 ) = MAXWRK
844 *
845 RETURN
846 *
847 9999 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
848 $ I6, ', JTYPE=', I6, ')' )
849 *
850 9998 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
851 $ I6, ', Input Example #', I2, ')' )
852 *
853 9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', I1, ' for eigenvalue ',
854 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
855 *
856 9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', I6, '.',
857 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
858 *
859 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
860 $ ' problem driver' )
861 *
862 9994 FORMAT( 'Input Example' )
863 *
864 9993 FORMAT( ' Matrix types: ', /
865 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
866 $ 'and B is the identity ', / ' matrix, ',
867 $ / ' 2: A and B are upper triangular matrices, ',
868 $ / ' 3: A and B are as type 2, but each second diagonal ',
869 $ 'block in A_11 and ', /
870 $ ' each third diaongal block in A_22 are 2x2 blocks,',
871 $ / ' 4: A and B are block diagonal matrices, ',
872 $ / ' 5: (A,B) has potentially close or common ',
873 $ 'eigenvalues.', / )
874 *
875 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
876 $ 'Q and Z are ', A, ',', / 19X,
877 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
878 $ / ' 1 = | A - Q S Z', A,
879 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
880 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
881 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
882 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
883 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
884 $ ' and diagonals of (S,T)', /
885 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
886 $ 'selected eigenvalues', /
887 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
888 $ 'DIFTRU/DIFEST > 10*THRESH',
889 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
890 $ 'when reordering fails', /
891 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
892 $ 'PLTRU/PLEST > THRESH', /
893 $ ' ( Test 10 is only for input examples )', / )
894 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
895 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
896 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
897 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.4 )
898 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
899 $ ' result ', I2, ' is', 0P, F8.2 )
900 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
901 $ ' result ', I2, ' is', 1P, E10.3 )
902 *
903 * End of SDRGSX
904 *
905 END
2 $ AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
3 $ WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
4 *
5 * -- LAPACK test routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
11 $ NOUT, NSIZE
12 REAL THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 INTEGER IWORK( * )
17 REAL A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
18 $ ALPHAR( * ), B( LDA, * ), BETA( * ),
19 $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
20 $ WORK( * ), Z( LDA, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem expert driver SGGESX.
28 *
29 * SGGESX factors A and B as Q S Z' and Q T Z', where ' means
30 * transpose, T is upper triangular, S is in generalized Schur form
31 * (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
32 * the 2x2 blocks corresponding to complex conjugate pairs of
33 * generalized eigenvalues), and Q and Z are orthogonal. It also
34 * computes the generalized eigenvalues (alpha(1),beta(1)), ...,
35 * (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
36 * characteristic equation
37 *
38 * det( A - w(j) B ) = 0
39 *
40 * Optionally it also reorders the eigenvalues so that a selected
41 * cluster of eigenvalues appears in the leading diagonal block of the
42 * Schur forms; computes a reciprocal condition number for the average
43 * of the selected eigenvalues; and computes a reciprocal condition
44 * number for the right and left deflating subspaces corresponding to
45 * the selected eigenvalues.
46 *
47 * When SDRGSX is called with NSIZE > 0, five (5) types of built-in
48 * matrix pairs are used to test the routine SGGESX.
49 *
50 * When SDRGSX is called with NSIZE = 0, it reads in test matrix data
51 * to test SGGESX.
52 *
53 * For each matrix pair, the following tests will be performed and
54 * compared with the threshhold THRESH except for the tests (7) and (9):
55 *
56 * (1) | A - Q S Z' | / ( |A| n ulp )
57 *
58 * (2) | B - Q T Z' | / ( |B| n ulp )
59 *
60 * (3) | I - QQ' | / ( n ulp )
61 *
62 * (4) | I - ZZ' | / ( n ulp )
63 *
64 * (5) if A is in Schur form (i.e. quasi-triangular form)
65 *
66 * (6) maximum over j of D(j) where:
67 *
68 * if alpha(j) is real:
69 * |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
70 * D(j) = ------------------------ + -----------------------
71 * max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
72 *
73 * if alpha(j) is complex:
74 * | det( s S - w T ) |
75 * D(j) = ---------------------------------------------------
76 * ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
77 *
78 * and S and T are here the 2 x 2 diagonal blocks of S and T
79 * corresponding to the j-th and j+1-th eigenvalues.
80 *
81 * (7) if sorting worked and SDIM is the number of eigenvalues
82 * which were selected.
83 *
84 * (8) the estimated value DIF does not differ from the true values of
85 * Difu and Difl more than a factor 10*THRESH. If the estimate DIF
86 * equals zero the corresponding true values of Difu and Difl
87 * should be less than EPS*norm(A, B). If the true value of Difu
88 * and Difl equal zero, the estimate DIF should be less than
89 * EPS*norm(A, B).
90 *
91 * (9) If INFO = N+3 is returned by SGGESX, the reordering "failed"
92 * and we check that DIF = PL = PR = 0 and that the true value of
93 * Difu and Difl is < EPS*norm(A, B). We count the events when
94 * INFO=N+3.
95 *
96 * For read-in test matrices, the above tests are run except that the
97 * exact value for DIF (and PL) is input data. Additionally, there is
98 * one more test run for read-in test matrices:
99 *
100 * (10) the estimated value PL does not differ from the true value of
101 * PLTRU more than a factor THRESH. If the estimate PL equals
102 * zero the corresponding true value of PLTRU should be less than
103 * EPS*norm(A, B). If the true value of PLTRU equal zero, the
104 * estimate PL should be less than EPS*norm(A, B).
105 *
106 * Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
107 * matrix pairs are generated and tested. NSIZE should be kept small.
108 *
109 * SVD (routine SGESVD) is used for computing the true value of DIF_u
110 * and DIF_l when testing the built-in test problems.
111 *
112 * Built-in Test Matrices
113 * ======================
114 *
115 * All built-in test matrices are the 2 by 2 block of triangular
116 * matrices
117 *
118 * A = [ A11 A12 ] and B = [ B11 B12 ]
119 * [ A22 ] [ B22 ]
120 *
121 * where for different type of A11 and A22 are given as the following.
122 * A12 and B12 are chosen so that the generalized Sylvester equation
123 *
124 * A11*R - L*A22 = -A12
125 * B11*R - L*B22 = -B12
126 *
127 * have prescribed solution R and L.
128 *
129 * Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
130 * B11 = I_m, B22 = I_k
131 * where J_k(a,b) is the k-by-k Jordan block with ``a'' on
132 * diagonal and ``b'' on superdiagonal.
133 *
134 * Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and
135 * B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
136 * A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
137 * B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
138 *
139 * Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each
140 * second diagonal block in A_11 and each third diagonal block
141 * in A_22 are made as 2 by 2 blocks.
142 *
143 * Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
144 * for i=1,...,m, j=1,...,m and
145 * A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
146 * for i=m+1,...,k, j=m+1,...,k
147 *
148 * Type 5: (A,B) and have potentially close or common eigenvalues and
149 * very large departure from block diagonality A_11 is chosen
150 * as the m x m leading submatrix of A_1:
151 * | 1 b |
152 * | -b 1 |
153 * | 1+d b |
154 * | -b 1+d |
155 * A_1 = | d 1 |
156 * | -1 d |
157 * | -d 1 |
158 * | -1 -d |
159 * | 1 |
160 * and A_22 is chosen as the k x k leading submatrix of A_2:
161 * | -1 b |
162 * | -b -1 |
163 * | 1-d b |
164 * | -b 1-d |
165 * A_2 = | d 1+b |
166 * | -1-b d |
167 * | -d 1+b |
168 * | -1+b -d |
169 * | 1-d |
170 * and matrix B are chosen as identity matrices (see SLATM5).
171 *
172 *
173 * Arguments
174 * =========
175 *
176 * NSIZE (input) INTEGER
177 * The maximum size of the matrices to use. NSIZE >= 0.
178 * If NSIZE = 0, no built-in tests matrices are used, but
179 * read-in test matrices are used to test SGGESX.
180 *
181 * NCMAX (input) INTEGER
182 * Maximum allowable NMAX for generating Kroneker matrix
183 * in call to SLAKF2
184 *
185 * THRESH (input) REAL
186 * A test will count as "failed" if the "error", computed as
187 * described above, exceeds THRESH. Note that the error
188 * is scaled to be O(1), so THRESH should be a reasonably
189 * small multiple of 1, e.g., 10 or 100. In particular,
190 * it should not depend on the precision (single vs. double)
191 * or the size of the matrix. THRESH >= 0.
192 *
193 * NIN (input) INTEGER
194 * The FORTRAN unit number for reading in the data file of
195 * problems to solve.
196 *
197 * NOUT (input) INTEGER
198 * The FORTRAN unit number for printing out error messages
199 * (e.g., if a routine returns IINFO not equal to 0.)
200 *
201 * A (workspace) REAL array, dimension (LDA, NSIZE)
202 * Used to store the matrix whose eigenvalues are to be
203 * computed. On exit, A contains the last matrix actually used.
204 *
205 * LDA (input) INTEGER
206 * The leading dimension of A, B, AI, BI, Z and Q,
207 * LDA >= max( 1, NSIZE ). For the read-in test,
208 * LDA >= max( 1, N ), N is the size of the test matrices.
209 *
210 * B (workspace) REAL array, dimension (LDA, NSIZE)
211 * Used to store the matrix whose eigenvalues are to be
212 * computed. On exit, B contains the last matrix actually used.
213 *
214 * AI (workspace) REAL array, dimension (LDA, NSIZE)
215 * Copy of A, modified by SGGESX.
216 *
217 * BI (workspace) REAL array, dimension (LDA, NSIZE)
218 * Copy of B, modified by SGGESX.
219 *
220 * Z (workspace) REAL array, dimension (LDA, NSIZE)
221 * Z holds the left Schur vectors computed by SGGESX.
222 *
223 * Q (workspace) REAL array, dimension (LDA, NSIZE)
224 * Q holds the right Schur vectors computed by SGGESX.
225 *
226 * ALPHAR (workspace) REAL array, dimension (NSIZE)
227 * ALPHAI (workspace) REAL array, dimension (NSIZE)
228 * BETA (workspace) REAL array, dimension (NSIZE)
229 * On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
230 *
231 * C (workspace) REAL array, dimension (LDC, LDC)
232 * Store the matrix generated by subroutine SLAKF2, this is the
233 * matrix formed by Kronecker products used for estimating
234 * DIF.
235 *
236 * LDC (input) INTEGER
237 * The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
238 *
239 * S (workspace) REAL array, dimension (LDC)
240 * Singular values of C
241 *
242 * WORK (workspace) REAL array, dimension (LWORK)
243 *
244 * LWORK (input) INTEGER
245 * The dimension of the array WORK.
246 * LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
247 *
248 * IWORK (workspace) INTEGER array, dimension (LIWORK)
249 *
250 * LIWORK (input) INTEGER
251 * The dimension of the array IWORK. LIWORK >= NSIZE + 6.
252 *
253 * BWORK (workspace) LOGICAL array, dimension (LDA)
254 *
255 * INFO (output) INTEGER
256 * = 0: successful exit
257 * < 0: if INFO = -i, the i-th argument had an illegal value.
258 * > 0: A routine returned an error code.
259 *
260 * =====================================================================
261 *
262 * .. Parameters ..
263 REAL ZERO, ONE, TEN
264 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
265 * ..
266 * .. Local Scalars ..
267 LOGICAL ILABAD
268 CHARACTER SENSE
269 INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
270 $ MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
271 $ PRTYPE, QBA, QBB
272 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
273 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
274 * ..
275 * .. Local Arrays ..
276 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
277 * ..
278 * .. External Functions ..
279 LOGICAL SLCTSX
280 INTEGER ILAENV
281 REAL SLAMCH, SLANGE
282 EXTERNAL SLCTSX, ILAENV, SLAMCH, SLANGE
283 * ..
284 * .. External Subroutines ..
285 EXTERNAL ALASVM, SGESVD, SGET51, SGET53, SGGESX, SLABAD,
286 $ SLACPY, SLAKF2, SLASET, SLATM5, XERBLA
287 * ..
288 * .. Intrinsic Functions ..
289 INTRINSIC ABS, MAX, SQRT
290 * ..
291 * .. Scalars in Common ..
292 LOGICAL FS
293 INTEGER K, M, MPLUSN, N
294 * ..
295 * .. Common blocks ..
296 COMMON / MN / M, N, MPLUSN, K, FS
297 * ..
298 * .. Executable Statements ..
299 *
300 * Check for errors
301 *
302 IF( NSIZE.LT.0 ) THEN
303 INFO = -1
304 ELSE IF( THRESH.LT.ZERO ) THEN
305 INFO = -2
306 ELSE IF( NIN.LE.0 ) THEN
307 INFO = -3
308 ELSE IF( NOUT.LE.0 ) THEN
309 INFO = -4
310 ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
311 INFO = -6
312 ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
313 INFO = -17
314 ELSE IF( LIWORK.LT.NSIZE+6 ) THEN
315 INFO = -21
316 END IF
317 *
318 * Compute workspace
319 * (Note: Comments in the code beginning "Workspace:" describe the
320 * minimal amount of workspace needed at that point in the code,
321 * as well as the preferred amount for good performance.
322 * NB refers to the optimal block size for the immediately
323 * following subroutine, as returned by ILAENV.)
324 *
325 MINWRK = 1
326 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
327 c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
328 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
329 *
330 * workspace for sggesx
331 *
332 MAXWRK = 9*( NSIZE+1 ) + NSIZE*
333 $ ILAENV( 1, 'SGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
334 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
335 $ ILAENV( 1, 'SORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
336 *
337 * workspace for sgesvd
338 *
339 BDSPAC = 5*NSIZE*NSIZE / 2
340 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
341 $ ILAENV( 1, 'SGEBRD', ' ', NSIZE*NSIZE / 2,
342 $ NSIZE*NSIZE / 2, -1, -1 ) )
343 MAXWRK = MAX( MAXWRK, BDSPAC )
344 *
345 MAXWRK = MAX( MAXWRK, MINWRK )
346 *
347 WORK( 1 ) = MAXWRK
348 END IF
349 *
350 IF( LWORK.LT.MINWRK )
351 $ INFO = -19
352 *
353 IF( INFO.NE.0 ) THEN
354 CALL XERBLA( 'SDRGSX', -INFO )
355 RETURN
356 END IF
357 *
358 * Important constants
359 *
360 ULP = SLAMCH( 'P' )
361 ULPINV = ONE / ULP
362 SMLNUM = SLAMCH( 'S' ) / ULP
363 BIGNUM = ONE / SMLNUM
364 CALL SLABAD( SMLNUM, BIGNUM )
365 THRSH2 = TEN*THRESH
366 NTESTT = 0
367 NERRS = 0
368 *
369 * Go to the tests for read-in matrix pairs
370 *
371 IFUNC = 0
372 IF( NSIZE.EQ.0 )
373 $ GO TO 70
374 *
375 * Test the built-in matrix pairs.
376 * Loop over different functions (IFUNC) of SGGESX, types (PRTYPE)
377 * of test matrices, different size (M+N)
378 *
379 PRTYPE = 0
380 QBA = 3
381 QBB = 4
382 WEIGHT = SQRT( ULP )
383 *
384 DO 60 IFUNC = 0, 3
385 DO 50 PRTYPE = 1, 5
386 DO 40 M = 1, NSIZE - 1
387 DO 30 N = 1, NSIZE - M
388 *
389 WEIGHT = ONE / WEIGHT
390 MPLUSN = M + N
391 *
392 * Generate test matrices
393 *
394 FS = .TRUE.
395 K = 0
396 *
397 CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
398 $ LDA )
399 CALL SLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
400 $ LDA )
401 *
402 CALL SLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
403 $ LDA, AI( 1, M+1 ), LDA, BI, LDA,
404 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
405 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
406 *
407 * Compute the Schur factorization and swapping the
408 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
409 * Swapping is accomplished via the function SLCTSX
410 * which is supplied below.
411 *
412 IF( IFUNC.EQ.0 ) THEN
413 SENSE = 'N'
414 ELSE IF( IFUNC.EQ.1 ) THEN
415 SENSE = 'E'
416 ELSE IF( IFUNC.EQ.2 ) THEN
417 SENSE = 'V'
418 ELSE IF( IFUNC.EQ.3 ) THEN
419 SENSE = 'B'
420 END IF
421 *
422 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
423 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
424 *
425 CALL SGGESX( 'V', 'V', 'S', SLCTSX, SENSE, MPLUSN, AI,
426 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
427 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
428 $ IWORK, LIWORK, BWORK, LINFO )
429 *
430 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
431 RESULT( 1 ) = ULPINV
432 WRITE( NOUT, FMT = 9999 )'SGGESX', LINFO, MPLUSN,
433 $ PRTYPE
434 INFO = LINFO
435 GO TO 30
436 END IF
437 *
438 * Compute the norm(A, B)
439 *
440 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
441 $ MPLUSN )
442 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
443 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
444 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
445 $ WORK )
446 *
447 * Do tests (1) to (4)
448 *
449 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
450 $ LDA, WORK, RESULT( 1 ) )
451 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
452 $ LDA, WORK, RESULT( 2 ) )
453 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
454 $ LDA, WORK, RESULT( 3 ) )
455 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
456 $ LDA, WORK, RESULT( 4 ) )
457 NTEST = 4
458 *
459 * Do tests (5) and (6): check Schur form of A and
460 * compare eigenvalues with diagonals.
461 *
462 TEMP1 = ZERO
463 RESULT( 5 ) = ZERO
464 RESULT( 6 ) = ZERO
465 *
466 DO 10 J = 1, MPLUSN
467 ILABAD = .FALSE.
468 IF( ALPHAI( J ).EQ.ZERO ) THEN
469 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
470 $ MAX( SMLNUM, ABS( ALPHAR( J ) ),
471 $ ABS( AI( J, J ) ) )+
472 $ ABS( BETA( J )-BI( J, J ) ) /
473 $ MAX( SMLNUM, ABS( BETA( J ) ),
474 $ ABS( BI( J, J ) ) ) ) / ULP
475 IF( J.LT.MPLUSN ) THEN
476 IF( AI( J+1, J ).NE.ZERO ) THEN
477 ILABAD = .TRUE.
478 RESULT( 5 ) = ULPINV
479 END IF
480 END IF
481 IF( J.GT.1 ) THEN
482 IF( AI( J, J-1 ).NE.ZERO ) THEN
483 ILABAD = .TRUE.
484 RESULT( 5 ) = ULPINV
485 END IF
486 END IF
487 ELSE
488 IF( ALPHAI( J ).GT.ZERO ) THEN
489 I1 = J
490 ELSE
491 I1 = J - 1
492 END IF
493 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
494 ILABAD = .TRUE.
495 ELSE IF( I1.LT.MPLUSN-1 ) THEN
496 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
497 ILABAD = .TRUE.
498 RESULT( 5 ) = ULPINV
499 END IF
500 ELSE IF( I1.GT.1 ) THEN
501 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
502 ILABAD = .TRUE.
503 RESULT( 5 ) = ULPINV
504 END IF
505 END IF
506 IF( .NOT.ILABAD ) THEN
507 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
508 $ LDA, BETA( J ), ALPHAR( J ),
509 $ ALPHAI( J ), TEMP2, IINFO )
510 IF( IINFO.GE.3 ) THEN
511 WRITE( NOUT, FMT = 9997 )IINFO, J,
512 $ MPLUSN, PRTYPE
513 INFO = ABS( IINFO )
514 END IF
515 ELSE
516 TEMP2 = ULPINV
517 END IF
518 END IF
519 TEMP1 = MAX( TEMP1, TEMP2 )
520 IF( ILABAD ) THEN
521 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
522 END IF
523 10 CONTINUE
524 RESULT( 6 ) = TEMP1
525 NTEST = NTEST + 2
526 *
527 * Test (7) (if sorting worked)
528 *
529 RESULT( 7 ) = ZERO
530 IF( LINFO.EQ.MPLUSN+3 ) THEN
531 RESULT( 7 ) = ULPINV
532 ELSE IF( MM.NE.N ) THEN
533 RESULT( 7 ) = ULPINV
534 END IF
535 NTEST = NTEST + 1
536 *
537 * Test (8): compare the estimated value DIF and its
538 * value. first, compute the exact DIF.
539 *
540 RESULT( 8 ) = ZERO
541 MN2 = MM*( MPLUSN-MM )*2
542 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
543 *
544 * Note: for either following two causes, there are
545 * almost same number of test cases fail the test.
546 *
547 CALL SLAKF2( MM, MPLUSN-MM, AI, LDA,
548 $ AI( MM+1, MM+1 ), BI,
549 $ BI( MM+1, MM+1 ), C, LDC )
550 *
551 CALL SGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
552 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
553 $ INFO )
554 DIFTRU = S( MN2 )
555 *
556 IF( DIFEST( 2 ).EQ.ZERO ) THEN
557 IF( DIFTRU.GT.ABNRM*ULP )
558 $ RESULT( 8 ) = ULPINV
559 ELSE IF( DIFTRU.EQ.ZERO ) THEN
560 IF( DIFEST( 2 ).GT.ABNRM*ULP )
561 $ RESULT( 8 ) = ULPINV
562 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
563 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
564 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
565 $ DIFEST( 2 ) / DIFTRU )
566 END IF
567 NTEST = NTEST + 1
568 END IF
569 *
570 * Test (9)
571 *
572 RESULT( 9 ) = ZERO
573 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
574 IF( DIFTRU.GT.ABNRM*ULP )
575 $ RESULT( 9 ) = ULPINV
576 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
577 $ RESULT( 9 ) = ULPINV
578 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
579 $ RESULT( 9 ) = ULPINV
580 NTEST = NTEST + 1
581 END IF
582 *
583 NTESTT = NTESTT + NTEST
584 *
585 * Print out tests which fail.
586 *
587 DO 20 J = 1, 9
588 IF( RESULT( J ).GE.THRESH ) THEN
589 *
590 * If this is the first test to fail,
591 * print a header to the data file.
592 *
593 IF( NERRS.EQ.0 ) THEN
594 WRITE( NOUT, FMT = 9995 )'SGX'
595 *
596 * Matrix types
597 *
598 WRITE( NOUT, FMT = 9993 )
599 *
600 * Tests performed
601 *
602 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
603 $ 'transpose', ( '''', I = 1, 4 )
604 *
605 END IF
606 NERRS = NERRS + 1
607 IF( RESULT( J ).LT.10000.0 ) THEN
608 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
609 $ WEIGHT, M, J, RESULT( J )
610 ELSE
611 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
612 $ WEIGHT, M, J, RESULT( J )
613 END IF
614 END IF
615 20 CONTINUE
616 *
617 30 CONTINUE
618 40 CONTINUE
619 50 CONTINUE
620 60 CONTINUE
621 *
622 GO TO 150
623 *
624 70 CONTINUE
625 *
626 * Read in data from file to check accuracy of condition estimation
627 * Read input data until N=0
628 *
629 NPTKNT = 0
630 *
631 80 CONTINUE
632 READ( NIN, FMT = *, END = 140 )MPLUSN
633 IF( MPLUSN.EQ.0 )
634 $ GO TO 140
635 READ( NIN, FMT = *, END = 140 )N
636 DO 90 I = 1, MPLUSN
637 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
638 90 CONTINUE
639 DO 100 I = 1, MPLUSN
640 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
641 100 CONTINUE
642 READ( NIN, FMT = * )PLTRU, DIFTRU
643 *
644 NPTKNT = NPTKNT + 1
645 FS = .TRUE.
646 K = 0
647 M = MPLUSN - N
648 *
649 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
650 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
651 *
652 * Compute the Schur factorization while swaping the
653 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
654 *
655 CALL SGGESX( 'V', 'V', 'S', SLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
656 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
657 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
658 *
659 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
660 RESULT( 1 ) = ULPINV
661 WRITE( NOUT, FMT = 9998 )'SGGESX', LINFO, MPLUSN, NPTKNT
662 GO TO 130
663 END IF
664 *
665 * Compute the norm(A, B)
666 * (should this be norm of (A,B) or (AI,BI)?)
667 *
668 CALL SLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
669 CALL SLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
670 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
671 ABNRM = SLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
672 *
673 * Do tests (1) to (4)
674 *
675 CALL SGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
676 $ RESULT( 1 ) )
677 CALL SGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
678 $ RESULT( 2 ) )
679 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
680 $ RESULT( 3 ) )
681 CALL SGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
682 $ RESULT( 4 ) )
683 *
684 * Do tests (5) and (6): check Schur form of A and compare
685 * eigenvalues with diagonals.
686 *
687 NTEST = 6
688 TEMP1 = ZERO
689 RESULT( 5 ) = ZERO
690 RESULT( 6 ) = ZERO
691 *
692 DO 110 J = 1, MPLUSN
693 ILABAD = .FALSE.
694 IF( ALPHAI( J ).EQ.ZERO ) THEN
695 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
696 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
697 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
698 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
699 $ / ULP
700 IF( J.LT.MPLUSN ) THEN
701 IF( AI( J+1, J ).NE.ZERO ) THEN
702 ILABAD = .TRUE.
703 RESULT( 5 ) = ULPINV
704 END IF
705 END IF
706 IF( J.GT.1 ) THEN
707 IF( AI( J, J-1 ).NE.ZERO ) THEN
708 ILABAD = .TRUE.
709 RESULT( 5 ) = ULPINV
710 END IF
711 END IF
712 ELSE
713 IF( ALPHAI( J ).GT.ZERO ) THEN
714 I1 = J
715 ELSE
716 I1 = J - 1
717 END IF
718 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
719 ILABAD = .TRUE.
720 ELSE IF( I1.LT.MPLUSN-1 ) THEN
721 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
722 ILABAD = .TRUE.
723 RESULT( 5 ) = ULPINV
724 END IF
725 ELSE IF( I1.GT.1 ) THEN
726 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
727 ILABAD = .TRUE.
728 RESULT( 5 ) = ULPINV
729 END IF
730 END IF
731 IF( .NOT.ILABAD ) THEN
732 CALL SGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
733 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
734 $ IINFO )
735 IF( IINFO.GE.3 ) THEN
736 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
737 INFO = ABS( IINFO )
738 END IF
739 ELSE
740 TEMP2 = ULPINV
741 END IF
742 END IF
743 TEMP1 = MAX( TEMP1, TEMP2 )
744 IF( ILABAD ) THEN
745 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
746 END IF
747 110 CONTINUE
748 RESULT( 6 ) = TEMP1
749 *
750 * Test (7) (if sorting worked) <--------- need to be checked.
751 *
752 NTEST = 7
753 RESULT( 7 ) = ZERO
754 IF( LINFO.EQ.MPLUSN+3 )
755 $ RESULT( 7 ) = ULPINV
756 *
757 * Test (8): compare the estimated value of DIF and its true value.
758 *
759 NTEST = 8
760 RESULT( 8 ) = ZERO
761 IF( DIFEST( 2 ).EQ.ZERO ) THEN
762 IF( DIFTRU.GT.ABNRM*ULP )
763 $ RESULT( 8 ) = ULPINV
764 ELSE IF( DIFTRU.EQ.ZERO ) THEN
765 IF( DIFEST( 2 ).GT.ABNRM*ULP )
766 $ RESULT( 8 ) = ULPINV
767 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
768 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
769 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
770 END IF
771 *
772 * Test (9)
773 *
774 NTEST = 9
775 RESULT( 9 ) = ZERO
776 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
777 IF( DIFTRU.GT.ABNRM*ULP )
778 $ RESULT( 9 ) = ULPINV
779 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
780 $ RESULT( 9 ) = ULPINV
781 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
782 $ RESULT( 9 ) = ULPINV
783 END IF
784 *
785 * Test (10): compare the estimated value of PL and it true value.
786 *
787 NTEST = 10
788 RESULT( 10 ) = ZERO
789 IF( PL( 1 ).EQ.ZERO ) THEN
790 IF( PLTRU.GT.ABNRM*ULP )
791 $ RESULT( 10 ) = ULPINV
792 ELSE IF( PLTRU.EQ.ZERO ) THEN
793 IF( PL( 1 ).GT.ABNRM*ULP )
794 $ RESULT( 10 ) = ULPINV
795 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
796 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
797 RESULT( 10 ) = ULPINV
798 END IF
799 *
800 NTESTT = NTESTT + NTEST
801 *
802 * Print out tests which fail.
803 *
804 DO 120 J = 1, NTEST
805 IF( RESULT( J ).GE.THRESH ) THEN
806 *
807 * If this is the first test to fail,
808 * print a header to the data file.
809 *
810 IF( NERRS.EQ.0 ) THEN
811 WRITE( NOUT, FMT = 9995 )'SGX'
812 *
813 * Matrix types
814 *
815 WRITE( NOUT, FMT = 9994 )
816 *
817 * Tests performed
818 *
819 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
820 $ 'transpose', ( '''', I = 1, 4 )
821 *
822 END IF
823 NERRS = NERRS + 1
824 IF( RESULT( J ).LT.10000.0 ) THEN
825 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
826 ELSE
827 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
828 END IF
829 END IF
830 *
831 120 CONTINUE
832 *
833 130 CONTINUE
834 GO TO 80
835 140 CONTINUE
836 *
837 150 CONTINUE
838 *
839 * Summary
840 *
841 CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
842 *
843 WORK( 1 ) = MAXWRK
844 *
845 RETURN
846 *
847 9999 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
848 $ I6, ', JTYPE=', I6, ')' )
849 *
850 9998 FORMAT( ' SDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
851 $ I6, ', Input Example #', I2, ')' )
852 *
853 9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', I1, ' for eigenvalue ',
854 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
855 *
856 9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', I6, '.',
857 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
858 *
859 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
860 $ ' problem driver' )
861 *
862 9994 FORMAT( 'Input Example' )
863 *
864 9993 FORMAT( ' Matrix types: ', /
865 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
866 $ 'and B is the identity ', / ' matrix, ',
867 $ / ' 2: A and B are upper triangular matrices, ',
868 $ / ' 3: A and B are as type 2, but each second diagonal ',
869 $ 'block in A_11 and ', /
870 $ ' each third diaongal block in A_22 are 2x2 blocks,',
871 $ / ' 4: A and B are block diagonal matrices, ',
872 $ / ' 5: (A,B) has potentially close or common ',
873 $ 'eigenvalues.', / )
874 *
875 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
876 $ 'Q and Z are ', A, ',', / 19X,
877 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
878 $ / ' 1 = | A - Q S Z', A,
879 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
880 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
881 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
882 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
883 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
884 $ ' and diagonals of (S,T)', /
885 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
886 $ 'selected eigenvalues', /
887 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
888 $ 'DIFTRU/DIFEST > 10*THRESH',
889 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
890 $ 'when reordering fails', /
891 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
892 $ 'PLTRU/PLEST > THRESH', /
893 $ ' ( Test 10 is only for input examples )', / )
894 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
895 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
896 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.4,
897 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.4 )
898 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
899 $ ' result ', I2, ' is', 0P, F8.2 )
900 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
901 $ ' result ', I2, ' is', 1P, E10.3 )
902 *
903 * End of SDRGSX
904 *
905 END