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