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