1 SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
2 $ 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 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 INTEGER IWORK( * )
17 DOUBLE PRECISION 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 * DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem expert driver DGGESX.
28 *
29 * DGGESX 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 DDRGSX is called with NSIZE > 0, five (5) types of built-in
48 * matrix pairs are used to test the routine DGGESX.
49 *
50 * When DDRGSX is called with NSIZE = 0, it reads in test matrix data
51 * to test DGGESX.
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 DGGESX, 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 DGESVD) 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 DLATM5).
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 DGGESX.
180 *
181 * NCMAX (input) INTEGER
182 * Maximum allowable NMAX for generating Kroneker matrix
183 * in call to DLAKF2
184 *
185 * THRESH (input) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDA, NSIZE)
215 * Copy of A, modified by DGGESX.
216 *
217 * BI (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
218 * Copy of B, modified by DGGESX.
219 *
220 * Z (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
221 * Z holds the left Schur vectors computed by DGGESX.
222 *
223 * Q (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
224 * Q holds the right Schur vectors computed by DGGESX.
225 *
226 * ALPHAR (workspace) DOUBLE PRECISION array, dimension (NSIZE)
227 * ALPHAI (workspace) DOUBLE PRECISION array, dimension (NSIZE)
228 * BETA (workspace) DOUBLE PRECISION array, dimension (NSIZE)
229 * On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
230 *
231 * C (workspace) DOUBLE PRECISION array, dimension (LDC, LDC)
232 * Store the matrix generated by subroutine DLAKF2, 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) DOUBLE PRECISION array, dimension (LDC)
240 * Singular values of C
241 *
242 * WORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE, TEN
264 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+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 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
273 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
274 * ..
275 * .. Local Arrays ..
276 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
277 * ..
278 * .. External Functions ..
279 LOGICAL DLCTSX
280 INTEGER ILAENV
281 DOUBLE PRECISION DLAMCH, DLANGE
282 EXTERNAL DLCTSX, ILAENV, DLAMCH, DLANGE
283 * ..
284 * .. External Subroutines ..
285 EXTERNAL ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD,
286 $ DLACPY, DLAKF2, DLASET, DLATM5, 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 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
328 *
329 * workspace for sggesx
330 *
331 MAXWRK = 9*( NSIZE+1 ) + NSIZE*
332 $ ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
333 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
334 $ ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
335 *
336 * workspace for dgesvd
337 *
338 BDSPAC = 5*NSIZE*NSIZE / 2
339 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
340 $ ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2,
341 $ NSIZE*NSIZE / 2, -1, -1 ) )
342 MAXWRK = MAX( MAXWRK, BDSPAC )
343 *
344 MAXWRK = MAX( MAXWRK, MINWRK )
345 *
346 WORK( 1 ) = MAXWRK
347 END IF
348 *
349 IF( LWORK.LT.MINWRK )
350 $ INFO = -19
351 *
352 IF( INFO.NE.0 ) THEN
353 CALL XERBLA( 'DDRGSX', -INFO )
354 RETURN
355 END IF
356 *
357 * Important constants
358 *
359 ULP = DLAMCH( 'P' )
360 ULPINV = ONE / ULP
361 SMLNUM = DLAMCH( 'S' ) / ULP
362 BIGNUM = ONE / SMLNUM
363 CALL DLABAD( SMLNUM, BIGNUM )
364 THRSH2 = TEN*THRESH
365 NTESTT = 0
366 NERRS = 0
367 *
368 * Go to the tests for read-in matrix pairs
369 *
370 IFUNC = 0
371 IF( NSIZE.EQ.0 )
372 $ GO TO 70
373 *
374 * Test the built-in matrix pairs.
375 * Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
376 * of test matrices, different size (M+N)
377 *
378 PRTYPE = 0
379 QBA = 3
380 QBB = 4
381 WEIGHT = SQRT( ULP )
382 *
383 DO 60 IFUNC = 0, 3
384 DO 50 PRTYPE = 1, 5
385 DO 40 M = 1, NSIZE - 1
386 DO 30 N = 1, NSIZE - M
387 *
388 WEIGHT = ONE / WEIGHT
389 MPLUSN = M + N
390 *
391 * Generate test matrices
392 *
393 FS = .TRUE.
394 K = 0
395 *
396 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
397 $ LDA )
398 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
399 $ LDA )
400 *
401 CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
402 $ LDA, AI( 1, M+1 ), LDA, BI, LDA,
403 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
404 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
405 *
406 * Compute the Schur factorization and swapping the
407 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
408 * Swapping is accomplished via the function DLCTSX
409 * which is supplied below.
410 *
411 IF( IFUNC.EQ.0 ) THEN
412 SENSE = 'N'
413 ELSE IF( IFUNC.EQ.1 ) THEN
414 SENSE = 'E'
415 ELSE IF( IFUNC.EQ.2 ) THEN
416 SENSE = 'V'
417 ELSE IF( IFUNC.EQ.3 ) THEN
418 SENSE = 'B'
419 END IF
420 *
421 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
422 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
423 *
424 CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI,
425 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
426 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
427 $ IWORK, LIWORK, BWORK, LINFO )
428 *
429 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
430 RESULT( 1 ) = ULPINV
431 WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN,
432 $ PRTYPE
433 INFO = LINFO
434 GO TO 30
435 END IF
436 *
437 * Compute the norm(A, B)
438 *
439 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
440 $ MPLUSN )
441 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
442 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
443 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
444 $ WORK )
445 *
446 * Do tests (1) to (4)
447 *
448 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
449 $ LDA, WORK, RESULT( 1 ) )
450 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
451 $ LDA, WORK, RESULT( 2 ) )
452 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
453 $ LDA, WORK, RESULT( 3 ) )
454 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
455 $ LDA, WORK, RESULT( 4 ) )
456 NTEST = 4
457 *
458 * Do tests (5) and (6): check Schur form of A and
459 * compare eigenvalues with diagonals.
460 *
461 TEMP1 = ZERO
462 RESULT( 5 ) = ZERO
463 RESULT( 6 ) = ZERO
464 *
465 DO 10 J = 1, MPLUSN
466 ILABAD = .FALSE.
467 IF( ALPHAI( J ).EQ.ZERO ) THEN
468 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
469 $ MAX( SMLNUM, ABS( ALPHAR( J ) ),
470 $ ABS( AI( J, J ) ) )+
471 $ ABS( BETA( J )-BI( J, J ) ) /
472 $ MAX( SMLNUM, ABS( BETA( J ) ),
473 $ ABS( BI( J, J ) ) ) ) / ULP
474 IF( J.LT.MPLUSN ) THEN
475 IF( AI( J+1, J ).NE.ZERO ) THEN
476 ILABAD = .TRUE.
477 RESULT( 5 ) = ULPINV
478 END IF
479 END IF
480 IF( J.GT.1 ) THEN
481 IF( AI( J, J-1 ).NE.ZERO ) THEN
482 ILABAD = .TRUE.
483 RESULT( 5 ) = ULPINV
484 END IF
485 END IF
486 ELSE
487 IF( ALPHAI( J ).GT.ZERO ) THEN
488 I1 = J
489 ELSE
490 I1 = J - 1
491 END IF
492 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
493 ILABAD = .TRUE.
494 ELSE IF( I1.LT.MPLUSN-1 ) THEN
495 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
496 ILABAD = .TRUE.
497 RESULT( 5 ) = ULPINV
498 END IF
499 ELSE IF( I1.GT.1 ) THEN
500 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
501 ILABAD = .TRUE.
502 RESULT( 5 ) = ULPINV
503 END IF
504 END IF
505 IF( .NOT.ILABAD ) THEN
506 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
507 $ LDA, BETA( J ), ALPHAR( J ),
508 $ ALPHAI( J ), TEMP2, IINFO )
509 IF( IINFO.GE.3 ) THEN
510 WRITE( NOUT, FMT = 9997 )IINFO, J,
511 $ MPLUSN, PRTYPE
512 INFO = ABS( IINFO )
513 END IF
514 ELSE
515 TEMP2 = ULPINV
516 END IF
517 END IF
518 TEMP1 = MAX( TEMP1, TEMP2 )
519 IF( ILABAD ) THEN
520 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
521 END IF
522 10 CONTINUE
523 RESULT( 6 ) = TEMP1
524 NTEST = NTEST + 2
525 *
526 * Test (7) (if sorting worked)
527 *
528 RESULT( 7 ) = ZERO
529 IF( LINFO.EQ.MPLUSN+3 ) THEN
530 RESULT( 7 ) = ULPINV
531 ELSE IF( MM.NE.N ) THEN
532 RESULT( 7 ) = ULPINV
533 END IF
534 NTEST = NTEST + 1
535 *
536 * Test (8): compare the estimated value DIF and its
537 * value. first, compute the exact DIF.
538 *
539 RESULT( 8 ) = ZERO
540 MN2 = MM*( MPLUSN-MM )*2
541 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
542 *
543 * Note: for either following two causes, there are
544 * almost same number of test cases fail the test.
545 *
546 CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
547 $ AI( MM+1, MM+1 ), BI,
548 $ BI( MM+1, MM+1 ), C, LDC )
549 *
550 CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
551 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
552 $ INFO )
553 DIFTRU = S( MN2 )
554 *
555 IF( DIFEST( 2 ).EQ.ZERO ) THEN
556 IF( DIFTRU.GT.ABNRM*ULP )
557 $ RESULT( 8 ) = ULPINV
558 ELSE IF( DIFTRU.EQ.ZERO ) THEN
559 IF( DIFEST( 2 ).GT.ABNRM*ULP )
560 $ RESULT( 8 ) = ULPINV
561 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
562 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
563 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
564 $ DIFEST( 2 ) / DIFTRU )
565 END IF
566 NTEST = NTEST + 1
567 END IF
568 *
569 * Test (9)
570 *
571 RESULT( 9 ) = ZERO
572 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
573 IF( DIFTRU.GT.ABNRM*ULP )
574 $ RESULT( 9 ) = ULPINV
575 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
576 $ RESULT( 9 ) = ULPINV
577 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
578 $ RESULT( 9 ) = ULPINV
579 NTEST = NTEST + 1
580 END IF
581 *
582 NTESTT = NTESTT + NTEST
583 *
584 * Print out tests which fail.
585 *
586 DO 20 J = 1, 9
587 IF( RESULT( J ).GE.THRESH ) THEN
588 *
589 * If this is the first test to fail,
590 * print a header to the data file.
591 *
592 IF( NERRS.EQ.0 ) THEN
593 WRITE( NOUT, FMT = 9995 )'SGX'
594 *
595 * Matrix types
596 *
597 WRITE( NOUT, FMT = 9993 )
598 *
599 * Tests performed
600 *
601 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
602 $ 'transpose', ( '''', I = 1, 4 )
603 *
604 END IF
605 NERRS = NERRS + 1
606 IF( RESULT( J ).LT.10000.0D0 ) THEN
607 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
608 $ WEIGHT, M, J, RESULT( J )
609 ELSE
610 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
611 $ WEIGHT, M, J, RESULT( J )
612 END IF
613 END IF
614 20 CONTINUE
615 *
616 30 CONTINUE
617 40 CONTINUE
618 50 CONTINUE
619 60 CONTINUE
620 *
621 GO TO 150
622 *
623 70 CONTINUE
624 *
625 * Read in data from file to check accuracy of condition estimation
626 * Read input data until N=0
627 *
628 NPTKNT = 0
629 *
630 80 CONTINUE
631 READ( NIN, FMT = *, END = 140 )MPLUSN
632 IF( MPLUSN.EQ.0 )
633 $ GO TO 140
634 READ( NIN, FMT = *, END = 140 )N
635 DO 90 I = 1, MPLUSN
636 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
637 90 CONTINUE
638 DO 100 I = 1, MPLUSN
639 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
640 100 CONTINUE
641 READ( NIN, FMT = * )PLTRU, DIFTRU
642 *
643 NPTKNT = NPTKNT + 1
644 FS = .TRUE.
645 K = 0
646 M = MPLUSN - N
647 *
648 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
649 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
650 *
651 * Compute the Schur factorization while swaping the
652 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
653 *
654 CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
655 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
656 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
657 *
658 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
659 RESULT( 1 ) = ULPINV
660 WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
661 GO TO 130
662 END IF
663 *
664 * Compute the norm(A, B)
665 * (should this be norm of (A,B) or (AI,BI)?)
666 *
667 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
668 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
669 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
670 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
671 *
672 * Do tests (1) to (4)
673 *
674 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
675 $ RESULT( 1 ) )
676 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
677 $ RESULT( 2 ) )
678 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
679 $ RESULT( 3 ) )
680 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
681 $ RESULT( 4 ) )
682 *
683 * Do tests (5) and (6): check Schur form of A and compare
684 * eigenvalues with diagonals.
685 *
686 NTEST = 6
687 TEMP1 = ZERO
688 RESULT( 5 ) = ZERO
689 RESULT( 6 ) = ZERO
690 *
691 DO 110 J = 1, MPLUSN
692 ILABAD = .FALSE.
693 IF( ALPHAI( J ).EQ.ZERO ) THEN
694 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
695 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
696 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
697 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
698 $ / ULP
699 IF( J.LT.MPLUSN ) THEN
700 IF( AI( J+1, J ).NE.ZERO ) THEN
701 ILABAD = .TRUE.
702 RESULT( 5 ) = ULPINV
703 END IF
704 END IF
705 IF( J.GT.1 ) THEN
706 IF( AI( J, J-1 ).NE.ZERO ) THEN
707 ILABAD = .TRUE.
708 RESULT( 5 ) = ULPINV
709 END IF
710 END IF
711 ELSE
712 IF( ALPHAI( J ).GT.ZERO ) THEN
713 I1 = J
714 ELSE
715 I1 = J - 1
716 END IF
717 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
718 ILABAD = .TRUE.
719 ELSE IF( I1.LT.MPLUSN-1 ) THEN
720 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
721 ILABAD = .TRUE.
722 RESULT( 5 ) = ULPINV
723 END IF
724 ELSE IF( I1.GT.1 ) THEN
725 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
726 ILABAD = .TRUE.
727 RESULT( 5 ) = ULPINV
728 END IF
729 END IF
730 IF( .NOT.ILABAD ) THEN
731 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
732 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
733 $ IINFO )
734 IF( IINFO.GE.3 ) THEN
735 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
736 INFO = ABS( IINFO )
737 END IF
738 ELSE
739 TEMP2 = ULPINV
740 END IF
741 END IF
742 TEMP1 = MAX( TEMP1, TEMP2 )
743 IF( ILABAD ) THEN
744 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
745 END IF
746 110 CONTINUE
747 RESULT( 6 ) = TEMP1
748 *
749 * Test (7) (if sorting worked) <--------- need to be checked.
750 *
751 NTEST = 7
752 RESULT( 7 ) = ZERO
753 IF( LINFO.EQ.MPLUSN+3 )
754 $ RESULT( 7 ) = ULPINV
755 *
756 * Test (8): compare the estimated value of DIF and its true value.
757 *
758 NTEST = 8
759 RESULT( 8 ) = ZERO
760 IF( DIFEST( 2 ).EQ.ZERO ) THEN
761 IF( DIFTRU.GT.ABNRM*ULP )
762 $ RESULT( 8 ) = ULPINV
763 ELSE IF( DIFTRU.EQ.ZERO ) THEN
764 IF( DIFEST( 2 ).GT.ABNRM*ULP )
765 $ RESULT( 8 ) = ULPINV
766 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
767 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
768 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
769 END IF
770 *
771 * Test (9)
772 *
773 NTEST = 9
774 RESULT( 9 ) = ZERO
775 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
776 IF( DIFTRU.GT.ABNRM*ULP )
777 $ RESULT( 9 ) = ULPINV
778 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
779 $ RESULT( 9 ) = ULPINV
780 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
781 $ RESULT( 9 ) = ULPINV
782 END IF
783 *
784 * Test (10): compare the estimated value of PL and it true value.
785 *
786 NTEST = 10
787 RESULT( 10 ) = ZERO
788 IF( PL( 1 ).EQ.ZERO ) THEN
789 IF( PLTRU.GT.ABNRM*ULP )
790 $ RESULT( 10 ) = ULPINV
791 ELSE IF( PLTRU.EQ.ZERO ) THEN
792 IF( PL( 1 ).GT.ABNRM*ULP )
793 $ RESULT( 10 ) = ULPINV
794 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
795 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
796 RESULT( 10 ) = ULPINV
797 END IF
798 *
799 NTESTT = NTESTT + NTEST
800 *
801 * Print out tests which fail.
802 *
803 DO 120 J = 1, NTEST
804 IF( RESULT( J ).GE.THRESH ) THEN
805 *
806 * If this is the first test to fail,
807 * print a header to the data file.
808 *
809 IF( NERRS.EQ.0 ) THEN
810 WRITE( NOUT, FMT = 9995 )'SGX'
811 *
812 * Matrix types
813 *
814 WRITE( NOUT, FMT = 9994 )
815 *
816 * Tests performed
817 *
818 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
819 $ 'transpose', ( '''', I = 1, 4 )
820 *
821 END IF
822 NERRS = NERRS + 1
823 IF( RESULT( J ).LT.10000.0D0 ) THEN
824 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
825 ELSE
826 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
827 END IF
828 END IF
829 *
830 120 CONTINUE
831 *
832 130 CONTINUE
833 GO TO 80
834 140 CONTINUE
835 *
836 150 CONTINUE
837 *
838 * Summary
839 *
840 CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
841 *
842 WORK( 1 ) = MAXWRK
843 *
844 RETURN
845 *
846 9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
847 $ I6, ', JTYPE=', I6, ')' )
848 *
849 9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
850 $ I6, ', Input Example #', I2, ')' )
851 *
852 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
853 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
854 *
855 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
856 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
857 *
858 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
859 $ ' problem driver' )
860 *
861 9994 FORMAT( 'Input Example' )
862 *
863 9993 FORMAT( ' Matrix types: ', /
864 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
865 $ 'and B is the identity ', / ' matrix, ',
866 $ / ' 2: A and B are upper triangular matrices, ',
867 $ / ' 3: A and B are as type 2, but each second diagonal ',
868 $ 'block in A_11 and ', /
869 $ ' each third diaongal block in A_22 are 2x2 blocks,',
870 $ / ' 4: A and B are block diagonal matrices, ',
871 $ / ' 5: (A,B) has potentially close or common ',
872 $ 'eigenvalues.', / )
873 *
874 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
875 $ 'Q and Z are ', A, ',', / 19X,
876 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
877 $ / ' 1 = | A - Q S Z', A,
878 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
879 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
880 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
881 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
882 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
883 $ ' and diagonals of (S,T)', /
884 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
885 $ 'selected eigenvalues', /
886 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
887 $ 'DIFTRU/DIFEST > 10*THRESH',
888 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
889 $ 'when reordering fails', /
890 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
891 $ 'PLTRU/PLEST > THRESH', /
892 $ ' ( Test 10 is only for input examples )', / )
893 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
894 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
895 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
896 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.4 )
897 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
898 $ ' result ', I2, ' is', 0P, F8.2 )
899 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
900 $ ' result ', I2, ' is', 1P, D10.3 )
901 *
902 * End of DDRGSX
903 *
904 END
2 $ 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 DOUBLE PRECISION THRESH
13 * ..
14 * .. Array Arguments ..
15 LOGICAL BWORK( * )
16 INTEGER IWORK( * )
17 DOUBLE PRECISION 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 * DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
27 * problem expert driver DGGESX.
28 *
29 * DGGESX 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 DDRGSX is called with NSIZE > 0, five (5) types of built-in
48 * matrix pairs are used to test the routine DGGESX.
49 *
50 * When DDRGSX is called with NSIZE = 0, it reads in test matrix data
51 * to test DGGESX.
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 DGGESX, 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 DGESVD) 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 DLATM5).
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 DGGESX.
180 *
181 * NCMAX (input) INTEGER
182 * Maximum allowable NMAX for generating Kroneker matrix
183 * in call to DLAKF2
184 *
185 * THRESH (input) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDA, NSIZE)
215 * Copy of A, modified by DGGESX.
216 *
217 * BI (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
218 * Copy of B, modified by DGGESX.
219 *
220 * Z (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
221 * Z holds the left Schur vectors computed by DGGESX.
222 *
223 * Q (workspace) DOUBLE PRECISION array, dimension (LDA, NSIZE)
224 * Q holds the right Schur vectors computed by DGGESX.
225 *
226 * ALPHAR (workspace) DOUBLE PRECISION array, dimension (NSIZE)
227 * ALPHAI (workspace) DOUBLE PRECISION array, dimension (NSIZE)
228 * BETA (workspace) DOUBLE PRECISION array, dimension (NSIZE)
229 * On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
230 *
231 * C (workspace) DOUBLE PRECISION array, dimension (LDC, LDC)
232 * Store the matrix generated by subroutine DLAKF2, 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) DOUBLE PRECISION array, dimension (LDC)
240 * Singular values of C
241 *
242 * WORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE, TEN
264 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+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 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
273 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
274 * ..
275 * .. Local Arrays ..
276 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
277 * ..
278 * .. External Functions ..
279 LOGICAL DLCTSX
280 INTEGER ILAENV
281 DOUBLE PRECISION DLAMCH, DLANGE
282 EXTERNAL DLCTSX, ILAENV, DLAMCH, DLANGE
283 * ..
284 * .. External Subroutines ..
285 EXTERNAL ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD,
286 $ DLACPY, DLAKF2, DLASET, DLATM5, 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 MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
328 *
329 * workspace for sggesx
330 *
331 MAXWRK = 9*( NSIZE+1 ) + NSIZE*
332 $ ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
333 MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
334 $ ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
335 *
336 * workspace for dgesvd
337 *
338 BDSPAC = 5*NSIZE*NSIZE / 2
339 MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
340 $ ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2,
341 $ NSIZE*NSIZE / 2, -1, -1 ) )
342 MAXWRK = MAX( MAXWRK, BDSPAC )
343 *
344 MAXWRK = MAX( MAXWRK, MINWRK )
345 *
346 WORK( 1 ) = MAXWRK
347 END IF
348 *
349 IF( LWORK.LT.MINWRK )
350 $ INFO = -19
351 *
352 IF( INFO.NE.0 ) THEN
353 CALL XERBLA( 'DDRGSX', -INFO )
354 RETURN
355 END IF
356 *
357 * Important constants
358 *
359 ULP = DLAMCH( 'P' )
360 ULPINV = ONE / ULP
361 SMLNUM = DLAMCH( 'S' ) / ULP
362 BIGNUM = ONE / SMLNUM
363 CALL DLABAD( SMLNUM, BIGNUM )
364 THRSH2 = TEN*THRESH
365 NTESTT = 0
366 NERRS = 0
367 *
368 * Go to the tests for read-in matrix pairs
369 *
370 IFUNC = 0
371 IF( NSIZE.EQ.0 )
372 $ GO TO 70
373 *
374 * Test the built-in matrix pairs.
375 * Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
376 * of test matrices, different size (M+N)
377 *
378 PRTYPE = 0
379 QBA = 3
380 QBB = 4
381 WEIGHT = SQRT( ULP )
382 *
383 DO 60 IFUNC = 0, 3
384 DO 50 PRTYPE = 1, 5
385 DO 40 M = 1, NSIZE - 1
386 DO 30 N = 1, NSIZE - M
387 *
388 WEIGHT = ONE / WEIGHT
389 MPLUSN = M + N
390 *
391 * Generate test matrices
392 *
393 FS = .TRUE.
394 K = 0
395 *
396 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
397 $ LDA )
398 CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
399 $ LDA )
400 *
401 CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
402 $ LDA, AI( 1, M+1 ), LDA, BI, LDA,
403 $ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
404 $ Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
405 *
406 * Compute the Schur factorization and swapping the
407 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
408 * Swapping is accomplished via the function DLCTSX
409 * which is supplied below.
410 *
411 IF( IFUNC.EQ.0 ) THEN
412 SENSE = 'N'
413 ELSE IF( IFUNC.EQ.1 ) THEN
414 SENSE = 'E'
415 ELSE IF( IFUNC.EQ.2 ) THEN
416 SENSE = 'V'
417 ELSE IF( IFUNC.EQ.3 ) THEN
418 SENSE = 'B'
419 END IF
420 *
421 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
422 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
423 *
424 CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI,
425 $ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
426 $ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
427 $ IWORK, LIWORK, BWORK, LINFO )
428 *
429 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
430 RESULT( 1 ) = ULPINV
431 WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN,
432 $ PRTYPE
433 INFO = LINFO
434 GO TO 30
435 END IF
436 *
437 * Compute the norm(A, B)
438 *
439 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
440 $ MPLUSN )
441 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
442 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
443 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
444 $ WORK )
445 *
446 * Do tests (1) to (4)
447 *
448 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
449 $ LDA, WORK, RESULT( 1 ) )
450 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
451 $ LDA, WORK, RESULT( 2 ) )
452 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
453 $ LDA, WORK, RESULT( 3 ) )
454 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
455 $ LDA, WORK, RESULT( 4 ) )
456 NTEST = 4
457 *
458 * Do tests (5) and (6): check Schur form of A and
459 * compare eigenvalues with diagonals.
460 *
461 TEMP1 = ZERO
462 RESULT( 5 ) = ZERO
463 RESULT( 6 ) = ZERO
464 *
465 DO 10 J = 1, MPLUSN
466 ILABAD = .FALSE.
467 IF( ALPHAI( J ).EQ.ZERO ) THEN
468 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
469 $ MAX( SMLNUM, ABS( ALPHAR( J ) ),
470 $ ABS( AI( J, J ) ) )+
471 $ ABS( BETA( J )-BI( J, J ) ) /
472 $ MAX( SMLNUM, ABS( BETA( J ) ),
473 $ ABS( BI( J, J ) ) ) ) / ULP
474 IF( J.LT.MPLUSN ) THEN
475 IF( AI( J+1, J ).NE.ZERO ) THEN
476 ILABAD = .TRUE.
477 RESULT( 5 ) = ULPINV
478 END IF
479 END IF
480 IF( J.GT.1 ) THEN
481 IF( AI( J, J-1 ).NE.ZERO ) THEN
482 ILABAD = .TRUE.
483 RESULT( 5 ) = ULPINV
484 END IF
485 END IF
486 ELSE
487 IF( ALPHAI( J ).GT.ZERO ) THEN
488 I1 = J
489 ELSE
490 I1 = J - 1
491 END IF
492 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
493 ILABAD = .TRUE.
494 ELSE IF( I1.LT.MPLUSN-1 ) THEN
495 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
496 ILABAD = .TRUE.
497 RESULT( 5 ) = ULPINV
498 END IF
499 ELSE IF( I1.GT.1 ) THEN
500 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
501 ILABAD = .TRUE.
502 RESULT( 5 ) = ULPINV
503 END IF
504 END IF
505 IF( .NOT.ILABAD ) THEN
506 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
507 $ LDA, BETA( J ), ALPHAR( J ),
508 $ ALPHAI( J ), TEMP2, IINFO )
509 IF( IINFO.GE.3 ) THEN
510 WRITE( NOUT, FMT = 9997 )IINFO, J,
511 $ MPLUSN, PRTYPE
512 INFO = ABS( IINFO )
513 END IF
514 ELSE
515 TEMP2 = ULPINV
516 END IF
517 END IF
518 TEMP1 = MAX( TEMP1, TEMP2 )
519 IF( ILABAD ) THEN
520 WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
521 END IF
522 10 CONTINUE
523 RESULT( 6 ) = TEMP1
524 NTEST = NTEST + 2
525 *
526 * Test (7) (if sorting worked)
527 *
528 RESULT( 7 ) = ZERO
529 IF( LINFO.EQ.MPLUSN+3 ) THEN
530 RESULT( 7 ) = ULPINV
531 ELSE IF( MM.NE.N ) THEN
532 RESULT( 7 ) = ULPINV
533 END IF
534 NTEST = NTEST + 1
535 *
536 * Test (8): compare the estimated value DIF and its
537 * value. first, compute the exact DIF.
538 *
539 RESULT( 8 ) = ZERO
540 MN2 = MM*( MPLUSN-MM )*2
541 IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
542 *
543 * Note: for either following two causes, there are
544 * almost same number of test cases fail the test.
545 *
546 CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
547 $ AI( MM+1, MM+1 ), BI,
548 $ BI( MM+1, MM+1 ), C, LDC )
549 *
550 CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
551 $ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
552 $ INFO )
553 DIFTRU = S( MN2 )
554 *
555 IF( DIFEST( 2 ).EQ.ZERO ) THEN
556 IF( DIFTRU.GT.ABNRM*ULP )
557 $ RESULT( 8 ) = ULPINV
558 ELSE IF( DIFTRU.EQ.ZERO ) THEN
559 IF( DIFEST( 2 ).GT.ABNRM*ULP )
560 $ RESULT( 8 ) = ULPINV
561 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
562 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
563 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
564 $ DIFEST( 2 ) / DIFTRU )
565 END IF
566 NTEST = NTEST + 1
567 END IF
568 *
569 * Test (9)
570 *
571 RESULT( 9 ) = ZERO
572 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
573 IF( DIFTRU.GT.ABNRM*ULP )
574 $ RESULT( 9 ) = ULPINV
575 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
576 $ RESULT( 9 ) = ULPINV
577 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
578 $ RESULT( 9 ) = ULPINV
579 NTEST = NTEST + 1
580 END IF
581 *
582 NTESTT = NTESTT + NTEST
583 *
584 * Print out tests which fail.
585 *
586 DO 20 J = 1, 9
587 IF( RESULT( J ).GE.THRESH ) THEN
588 *
589 * If this is the first test to fail,
590 * print a header to the data file.
591 *
592 IF( NERRS.EQ.0 ) THEN
593 WRITE( NOUT, FMT = 9995 )'SGX'
594 *
595 * Matrix types
596 *
597 WRITE( NOUT, FMT = 9993 )
598 *
599 * Tests performed
600 *
601 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
602 $ 'transpose', ( '''', I = 1, 4 )
603 *
604 END IF
605 NERRS = NERRS + 1
606 IF( RESULT( J ).LT.10000.0D0 ) THEN
607 WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
608 $ WEIGHT, M, J, RESULT( J )
609 ELSE
610 WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
611 $ WEIGHT, M, J, RESULT( J )
612 END IF
613 END IF
614 20 CONTINUE
615 *
616 30 CONTINUE
617 40 CONTINUE
618 50 CONTINUE
619 60 CONTINUE
620 *
621 GO TO 150
622 *
623 70 CONTINUE
624 *
625 * Read in data from file to check accuracy of condition estimation
626 * Read input data until N=0
627 *
628 NPTKNT = 0
629 *
630 80 CONTINUE
631 READ( NIN, FMT = *, END = 140 )MPLUSN
632 IF( MPLUSN.EQ.0 )
633 $ GO TO 140
634 READ( NIN, FMT = *, END = 140 )N
635 DO 90 I = 1, MPLUSN
636 READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
637 90 CONTINUE
638 DO 100 I = 1, MPLUSN
639 READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
640 100 CONTINUE
641 READ( NIN, FMT = * )PLTRU, DIFTRU
642 *
643 NPTKNT = NPTKNT + 1
644 FS = .TRUE.
645 K = 0
646 M = MPLUSN - N
647 *
648 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
649 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
650 *
651 * Compute the Schur factorization while swaping the
652 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
653 *
654 CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
655 $ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
656 $ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
657 *
658 IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
659 RESULT( 1 ) = ULPINV
660 WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
661 GO TO 130
662 END IF
663 *
664 * Compute the norm(A, B)
665 * (should this be norm of (A,B) or (AI,BI)?)
666 *
667 CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
668 CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
669 $ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
670 ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
671 *
672 * Do tests (1) to (4)
673 *
674 CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
675 $ RESULT( 1 ) )
676 CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
677 $ RESULT( 2 ) )
678 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
679 $ RESULT( 3 ) )
680 CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
681 $ RESULT( 4 ) )
682 *
683 * Do tests (5) and (6): check Schur form of A and compare
684 * eigenvalues with diagonals.
685 *
686 NTEST = 6
687 TEMP1 = ZERO
688 RESULT( 5 ) = ZERO
689 RESULT( 6 ) = ZERO
690 *
691 DO 110 J = 1, MPLUSN
692 ILABAD = .FALSE.
693 IF( ALPHAI( J ).EQ.ZERO ) THEN
694 TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
695 $ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
696 $ J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
697 $ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
698 $ / ULP
699 IF( J.LT.MPLUSN ) THEN
700 IF( AI( J+1, J ).NE.ZERO ) THEN
701 ILABAD = .TRUE.
702 RESULT( 5 ) = ULPINV
703 END IF
704 END IF
705 IF( J.GT.1 ) THEN
706 IF( AI( J, J-1 ).NE.ZERO ) THEN
707 ILABAD = .TRUE.
708 RESULT( 5 ) = ULPINV
709 END IF
710 END IF
711 ELSE
712 IF( ALPHAI( J ).GT.ZERO ) THEN
713 I1 = J
714 ELSE
715 I1 = J - 1
716 END IF
717 IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
718 ILABAD = .TRUE.
719 ELSE IF( I1.LT.MPLUSN-1 ) THEN
720 IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
721 ILABAD = .TRUE.
722 RESULT( 5 ) = ULPINV
723 END IF
724 ELSE IF( I1.GT.1 ) THEN
725 IF( AI( I1, I1-1 ).NE.ZERO ) THEN
726 ILABAD = .TRUE.
727 RESULT( 5 ) = ULPINV
728 END IF
729 END IF
730 IF( .NOT.ILABAD ) THEN
731 CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
732 $ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
733 $ IINFO )
734 IF( IINFO.GE.3 ) THEN
735 WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
736 INFO = ABS( IINFO )
737 END IF
738 ELSE
739 TEMP2 = ULPINV
740 END IF
741 END IF
742 TEMP1 = MAX( TEMP1, TEMP2 )
743 IF( ILABAD ) THEN
744 WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
745 END IF
746 110 CONTINUE
747 RESULT( 6 ) = TEMP1
748 *
749 * Test (7) (if sorting worked) <--------- need to be checked.
750 *
751 NTEST = 7
752 RESULT( 7 ) = ZERO
753 IF( LINFO.EQ.MPLUSN+3 )
754 $ RESULT( 7 ) = ULPINV
755 *
756 * Test (8): compare the estimated value of DIF and its true value.
757 *
758 NTEST = 8
759 RESULT( 8 ) = ZERO
760 IF( DIFEST( 2 ).EQ.ZERO ) THEN
761 IF( DIFTRU.GT.ABNRM*ULP )
762 $ RESULT( 8 ) = ULPINV
763 ELSE IF( DIFTRU.EQ.ZERO ) THEN
764 IF( DIFEST( 2 ).GT.ABNRM*ULP )
765 $ RESULT( 8 ) = ULPINV
766 ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
767 $ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
768 RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
769 END IF
770 *
771 * Test (9)
772 *
773 NTEST = 9
774 RESULT( 9 ) = ZERO
775 IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
776 IF( DIFTRU.GT.ABNRM*ULP )
777 $ RESULT( 9 ) = ULPINV
778 IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
779 $ RESULT( 9 ) = ULPINV
780 IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
781 $ RESULT( 9 ) = ULPINV
782 END IF
783 *
784 * Test (10): compare the estimated value of PL and it true value.
785 *
786 NTEST = 10
787 RESULT( 10 ) = ZERO
788 IF( PL( 1 ).EQ.ZERO ) THEN
789 IF( PLTRU.GT.ABNRM*ULP )
790 $ RESULT( 10 ) = ULPINV
791 ELSE IF( PLTRU.EQ.ZERO ) THEN
792 IF( PL( 1 ).GT.ABNRM*ULP )
793 $ RESULT( 10 ) = ULPINV
794 ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
795 $ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
796 RESULT( 10 ) = ULPINV
797 END IF
798 *
799 NTESTT = NTESTT + NTEST
800 *
801 * Print out tests which fail.
802 *
803 DO 120 J = 1, NTEST
804 IF( RESULT( J ).GE.THRESH ) THEN
805 *
806 * If this is the first test to fail,
807 * print a header to the data file.
808 *
809 IF( NERRS.EQ.0 ) THEN
810 WRITE( NOUT, FMT = 9995 )'SGX'
811 *
812 * Matrix types
813 *
814 WRITE( NOUT, FMT = 9994 )
815 *
816 * Tests performed
817 *
818 WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
819 $ 'transpose', ( '''', I = 1, 4 )
820 *
821 END IF
822 NERRS = NERRS + 1
823 IF( RESULT( J ).LT.10000.0D0 ) THEN
824 WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
825 ELSE
826 WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
827 END IF
828 END IF
829 *
830 120 CONTINUE
831 *
832 130 CONTINUE
833 GO TO 80
834 140 CONTINUE
835 *
836 150 CONTINUE
837 *
838 * Summary
839 *
840 CALL ALASVM( 'SGX', NOUT, NERRS, NTESTT, 0 )
841 *
842 WORK( 1 ) = MAXWRK
843 *
844 RETURN
845 *
846 9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
847 $ I6, ', JTYPE=', I6, ')' )
848 *
849 9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
850 $ I6, ', Input Example #', I2, ')' )
851 *
852 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
853 $ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
854 *
855 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
856 $ / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
857 *
858 9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
859 $ ' problem driver' )
860 *
861 9994 FORMAT( 'Input Example' )
862 *
863 9993 FORMAT( ' Matrix types: ', /
864 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
865 $ 'and B is the identity ', / ' matrix, ',
866 $ / ' 2: A and B are upper triangular matrices, ',
867 $ / ' 3: A and B are as type 2, but each second diagonal ',
868 $ 'block in A_11 and ', /
869 $ ' each third diaongal block in A_22 are 2x2 blocks,',
870 $ / ' 4: A and B are block diagonal matrices, ',
871 $ / ' 5: (A,B) has potentially close or common ',
872 $ 'eigenvalues.', / )
873 *
874 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
875 $ 'Q and Z are ', A, ',', / 19X,
876 $ ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
877 $ / ' 1 = | A - Q S Z', A,
878 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
879 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
880 $ ' | / ( n ulp ) 4 = | I - ZZ', A,
881 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
882 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
883 $ ' and diagonals of (S,T)', /
884 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
885 $ 'selected eigenvalues', /
886 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
887 $ 'DIFTRU/DIFEST > 10*THRESH',
888 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
889 $ 'when reordering fails', /
890 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
891 $ 'PLTRU/PLEST > THRESH', /
892 $ ' ( Test 10 is only for input examples )', / )
893 9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
894 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
895 9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.4,
896 $ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.4 )
897 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
898 $ ' result ', I2, ' is', 0P, F8.2 )
899 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
900 $ ' result ', I2, ' is', 1P, D10.3 )
901 *
902 * End of DDRGSX
903 *
904 END