LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
cdrgsx.f
Go to the documentation of this file.
1 *> \brief \b CDRGSX
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 CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
12 * AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK,
13 * LWORK, RWORK, IWORK, LIWORK, BWORK, INFO )
14 *
15 * .. Scalar Arguments ..
16 * INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
17 * $ NOUT, NSIZE
18 * REAL THRESH
19 * ..
20 * .. Array Arguments ..
21 * LOGICAL BWORK( * )
22 * INTEGER IWORK( * )
23 * REAL RWORK( * ), S( * )
24 * COMPLEX 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 *> CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
37 *> problem expert driver CGGESX.
38 *>
39 *> CGGES 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 CDRGSX is called with NSIZE > 0, five (5) types of built-in
55 *> matrix pairs are used to test the routine CGGESX.
56 *>
57 *> When CDRGSX is called with NSIZE = 0, it reads in test matrix data
58 *> to test CGGESX.
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 CGGESX, 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 CGESVD) 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 SLATM5).
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 SGGESX.
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 CLAKF2
189 *> \endverbatim
190 *>
191 *> \param[in] THRESH
192 *> \verbatim
193 *> THRESH is REAL
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 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 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 array, dimension (LDA, NSIZE)
241 *> Copy of A, modified by CGGESX.
242 *> \endverbatim
243 *>
244 *> \param[out] BI
245 *> \verbatim
246 *> BI is COMPLEX array, dimension (LDA, NSIZE)
247 *> Copy of B, modified by CGGESX.
248 *> \endverbatim
249 *>
250 *> \param[out] Z
251 *> \verbatim
252 *> Z is COMPLEX array, dimension (LDA, NSIZE)
253 *> Z holds the left Schur vectors computed by CGGESX.
254 *> \endverbatim
255 *>
256 *> \param[out] Q
257 *> \verbatim
258 *> Q is COMPLEX array, dimension (LDA, NSIZE)
259 *> Q holds the right Schur vectors computed by CGGESX.
260 *> \endverbatim
261 *>
262 *> \param[out] ALPHA
263 *> \verbatim
264 *> ALPHA is COMPLEX array, dimension (NSIZE)
265 *> \endverbatim
266 *>
267 *> \param[out] BETA
268 *> \verbatim
269 *> BETA is COMPLEX 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 array, dimension (LDC, LDC)
277 *> Store the matrix generated by subroutine CLAKF2, 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 REAL array, dimension (LDC)
291 *> Singular values of C
292 *> \endverbatim
293 *>
294 *> \param[out] WORK
295 *> \verbatim
296 *> WORK is COMPLEX 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 REAL 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 complex_eig
346 *
347 * =====================================================================
348  SUBROUTINE cdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
349  $ ai, bi, z, q, alpha, beta, c, ldc, s, work,
350  $ lwork, 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  REAL thresh
361 * ..
362 * .. Array Arguments ..
363  LOGICAL bwork( * )
364  INTEGER iwork( * )
365  REAL rwork( * ), s( * )
366  COMPLEX 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  REAL zero, one, ten
376  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1 )
377  COMPLEX czero
378  parameter( czero = ( 0.0e+0, 0.0e+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  REAL abnrm, bignum, diftru, pltru, smlnum, temp1,
387  $ temp2, thrsh2, ulp, ulpinv, weight
388  COMPLEX x
389 * ..
390 * .. Local Arrays ..
391  REAL difest( 2 ), pl( 2 ), result( 10 )
392 * ..
393 * .. External Functions ..
394  LOGICAL clctsx
395  INTEGER ilaenv
396  REAL clange, slamch
397  EXTERNAL clctsx, ilaenv, clange, slamch
398 * ..
399 * .. External Subroutines ..
400  EXTERNAL alasvm, cgesvd, cget51, cggesx, clacpy, clakf2,
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, aimag, max, REAL, sqrt
412 * ..
413 * .. Statement Functions ..
414  REAL abs1
415 * ..
416 * .. Statement Function definitions ..
417  abs1( x ) = abs( REAL( X ) ) + abs( aimag( x ) )
418 * ..
419 * .. Executable Statements ..
420 *
421 * Check for errors
422 *
423  IF( nsize.LT.0 ) THEN
424  info = -1
425  ELSE IF( thresh.LT.zero ) THEN
426  info = -2
427  ELSE IF( nin.LE.0 ) THEN
428  info = -3
429  ELSE IF( nout.LE.0 ) THEN
430  info = -4
431  ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
432  info = -6
433  ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
434  info = -15
435  ELSE IF( liwork.LT.nsize+2 ) THEN
436  info = -21
437  END IF
438 *
439 * Compute workspace
440 * (Note: Comments in the code beginning "Workspace:" describe the
441 * minimal amount of workspace needed at that point in the code,
442 * as well as the preferred amount for good performance.
443 * NB refers to the optimal block size for the immediately
444 * following subroutine, as returned by ILAENV.)
445 *
446  minwrk = 1
447  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
448  minwrk = 3*nsize*nsize / 2
449 *
450 * workspace for cggesx
451 *
452  maxwrk = nsize*( 1+ilaenv( 1, 'CGEQRF', ' ', nsize, 1, nsize,
453  $ 0 ) )
454  maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1, 'CUNGQR', ' ',
455  $ nsize, 1, nsize, -1 ) ) )
456 *
457 * workspace for cgesvd
458 *
459  bdspac = 3*nsize*nsize / 2
460  maxwrk = max( maxwrk, nsize*nsize*
461  $ ( 1+ilaenv( 1, 'CGEBRD', ' ', nsize*nsize / 2,
462  $ nsize*nsize / 2, -1, -1 ) ) )
463  maxwrk = max( maxwrk, bdspac )
464 *
465  maxwrk = max( maxwrk, minwrk )
466 *
467  work( 1 ) = maxwrk
468  END IF
469 *
470  IF( lwork.LT.minwrk )
471  $ info = -18
472 *
473  IF( info.NE.0 ) THEN
474  CALL xerbla( 'CDRGSX', -info )
475  RETURN
476  END IF
477 *
478 * Important constants
479 *
480  ulp = slamch( 'P' )
481  ulpinv = one / ulp
482  smlnum = slamch( 'S' ) / ulp
483  bignum = one / smlnum
484  CALL slabad( smlnum, bignum )
485  thrsh2 = ten*thresh
486  ntestt = 0
487  nerrs = 0
488 *
489 * Go to the tests for read-in matrix pairs
490 *
491  ifunc = 0
492  IF( nsize.EQ.0 )
493  $ go to 70
494 *
495 * Test the built-in matrix pairs.
496 * Loop over different functions (IFUNC) of CGGESX, types (PRTYPE)
497 * of test matrices, different size (M+N)
498 *
499  prtype = 0
500  qba = 3
501  qbb = 4
502  weight = sqrt( ulp )
503 *
504  DO 60 ifunc = 0, 3
505  DO 50 prtype = 1, 5
506  DO 40 m = 1, nsize - 1
507  DO 30 n = 1, nsize - m
508 *
509  weight = one / weight
510  mplusn = m + n
511 *
512 * Generate test matrices
513 *
514  fs = .true.
515  k = 0
516 *
517  CALL claset( 'Full', mplusn, mplusn, czero, czero, ai,
518  $ lda )
519  CALL claset( 'Full', mplusn, mplusn, czero, czero, bi,
520  $ lda )
521 *
522  CALL clatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
523  $ lda, ai( 1, m+1 ), lda, bi, lda,
524  $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
525  $ q, lda, z, lda, weight, qba, qbb )
526 *
527 * Compute the Schur factorization and swapping the
528 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
529 * Swapping is accomplished via the function CLCTSX
530 * which is supplied below.
531 *
532  IF( ifunc.EQ.0 ) THEN
533  sense = 'N'
534  ELSE IF( ifunc.EQ.1 ) THEN
535  sense = 'E'
536  ELSE IF( ifunc.EQ.2 ) THEN
537  sense = 'V'
538  ELSE IF( ifunc.EQ.3 ) THEN
539  sense = 'B'
540  END IF
541 *
542  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
543  CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
544 *
545  CALL cggesx( 'V', 'V', 'S', clctsx, sense, mplusn, ai,
546  $ lda, bi, lda, mm, alpha, beta, q, lda, z,
547  $ lda, pl, difest, work, lwork, rwork,
548  $ iwork, liwork, bwork, linfo )
549 *
550  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
551  result( 1 ) = ulpinv
552  WRITE( nout, fmt = 9999 )'CGGESX', linfo, mplusn,
553  $ prtype
554  info = linfo
555  go to 30
556  END IF
557 *
558 * Compute the norm(A, B)
559 *
560  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work,
561  $ mplusn )
562  CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
563  $ work( mplusn*mplusn+1 ), mplusn )
564  abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn,
565  $ rwork )
566 *
567 * Do tests (1) to (4)
568 *
569  result( 2 ) = zero
570  CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
571  $ lda, work, rwork, result( 1 ) )
572  CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
573  $ lda, work, rwork, result( 2 ) )
574  CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
575  $ lda, work, rwork, result( 3 ) )
576  CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
577  $ lda, work, rwork, result( 4 ) )
578  ntest = 4
579 *
580 * Do tests (5) and (6): check Schur form of A and
581 * compare eigenvalues with diagonals.
582 *
583  temp1 = zero
584  result( 5 ) = zero
585  result( 6 ) = zero
586 *
587  DO 10 j = 1, mplusn
588  ilabad = .false.
589  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
590  $ max( smlnum, abs1( alpha( j ) ),
591  $ abs1( ai( j, j ) ) )+
592  $ abs1( beta( j )-bi( j, j ) ) /
593  $ max( smlnum, abs1( beta( j ) ),
594  $ abs1( bi( j, j ) ) ) ) / ulp
595  IF( j.LT.mplusn ) THEN
596  IF( ai( j+1, j ).NE.zero ) THEN
597  ilabad = .true.
598  result( 5 ) = ulpinv
599  END IF
600  END IF
601  IF( j.GT.1 ) THEN
602  IF( ai( j, j-1 ).NE.zero ) THEN
603  ilabad = .true.
604  result( 5 ) = ulpinv
605  END IF
606  END IF
607  temp1 = max( temp1, temp2 )
608  IF( ilabad ) THEN
609  WRITE( nout, fmt = 9997 )j, mplusn, prtype
610  END IF
611  10 CONTINUE
612  result( 6 ) = temp1
613  ntest = ntest + 2
614 *
615 * Test (7) (if sorting worked)
616 *
617  result( 7 ) = zero
618  IF( linfo.EQ.mplusn+3 ) THEN
619  result( 7 ) = ulpinv
620  ELSE IF( mm.NE.n ) THEN
621  result( 7 ) = ulpinv
622  END IF
623  ntest = ntest + 1
624 *
625 * Test (8): compare the estimated value DIF and its
626 * value. first, compute the exact DIF.
627 *
628  result( 8 ) = zero
629  mn2 = mm*( mplusn-mm )*2
630  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
631 *
632 * Note: for either following two cases, there are
633 * almost same number of test cases fail the test.
634 *
635  CALL clakf2( mm, mplusn-mm, ai, lda,
636  $ ai( mm+1, mm+1 ), bi,
637  $ bi( mm+1, mm+1 ), c, ldc )
638 *
639  CALL cgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
640  $ 1, work( 2 ), 1, work( 3 ), lwork-2,
641  $ rwork, info )
642  diftru = s( mn2 )
643 *
644  IF( difest( 2 ).EQ.zero ) THEN
645  IF( diftru.GT.abnrm*ulp )
646  $ result( 8 ) = ulpinv
647  ELSE IF( diftru.EQ.zero ) THEN
648  IF( difest( 2 ).GT.abnrm*ulp )
649  $ result( 8 ) = ulpinv
650  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
651  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
652  result( 8 ) = max( diftru / difest( 2 ),
653  $ difest( 2 ) / diftru )
654  END IF
655  ntest = ntest + 1
656  END IF
657 *
658 * Test (9)
659 *
660  result( 9 ) = zero
661  IF( linfo.EQ.( mplusn+2 ) ) THEN
662  IF( diftru.GT.abnrm*ulp )
663  $ result( 9 ) = ulpinv
664  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
665  $ result( 9 ) = ulpinv
666  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
667  $ result( 9 ) = ulpinv
668  ntest = ntest + 1
669  END IF
670 *
671  ntestt = ntestt + ntest
672 *
673 * Print out tests which fail.
674 *
675  DO 20 j = 1, 9
676  IF( result( j ).GE.thresh ) THEN
677 *
678 * If this is the first test to fail,
679 * print a header to the data file.
680 *
681  IF( nerrs.EQ.0 ) THEN
682  WRITE( nout, fmt = 9996 )'CGX'
683 *
684 * Matrix types
685 *
686  WRITE( nout, fmt = 9994 )
687 *
688 * Tests performed
689 *
690  WRITE( nout, fmt = 9993 )'unitary', '''',
691  $ 'transpose', ( '''', i = 1, 4 )
692 *
693  END IF
694  nerrs = nerrs + 1
695  IF( result( j ).LT.10000.0 ) THEN
696  WRITE( nout, fmt = 9992 )mplusn, prtype,
697  $ weight, m, j, result( j )
698  ELSE
699  WRITE( nout, fmt = 9991 )mplusn, prtype,
700  $ weight, m, j, result( j )
701  END IF
702  END IF
703  20 CONTINUE
704 *
705  30 CONTINUE
706  40 CONTINUE
707  50 CONTINUE
708  60 CONTINUE
709 *
710  go to 150
711 *
712  70 CONTINUE
713 *
714 * Read in data from file to check accuracy of condition estimation
715 * Read input data until N=0
716 *
717  nptknt = 0
718 *
719  80 CONTINUE
720  READ( nin, fmt = *, END = 140 )mplusn
721  IF( mplusn.EQ.0 )
722  $ go to 140
723  READ( nin, fmt = *, END = 140 )n
724  DO 90 i = 1, mplusn
725  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
726  90 CONTINUE
727  DO 100 i = 1, mplusn
728  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
729  100 CONTINUE
730  READ( nin, fmt = * )pltru, diftru
731 *
732  nptknt = nptknt + 1
733  fs = .true.
734  k = 0
735  m = mplusn - n
736 *
737  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
738  CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
739 *
740 * Compute the Schur factorization while swaping the
741 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
742 *
743  CALL cggesx( 'V', 'V', 'S', clctsx, 'B', mplusn, ai, lda, bi, lda,
744  $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
745  $ lwork, rwork, iwork, liwork, bwork, linfo )
746 *
747  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
748  result( 1 ) = ulpinv
749  WRITE( nout, fmt = 9998 )'CGGESX', linfo, mplusn, nptknt
750  go to 130
751  END IF
752 *
753 * Compute the norm(A, B)
754 * (should this be norm of (A,B) or (AI,BI)?)
755 *
756  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
757  CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
758  $ work( mplusn*mplusn+1 ), mplusn )
759  abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
760 *
761 * Do tests (1) to (4)
762 *
763  CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
764  $ rwork, result( 1 ) )
765  CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
766  $ rwork, result( 2 ) )
767  CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
768  $ rwork, result( 3 ) )
769  CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
770  $ rwork, result( 4 ) )
771 *
772 * Do tests (5) and (6): check Schur form of A and compare
773 * eigenvalues with diagonals.
774 *
775  ntest = 6
776  temp1 = zero
777  result( 5 ) = zero
778  result( 6 ) = zero
779 *
780  DO 110 j = 1, mplusn
781  ilabad = .false.
782  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
783  $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
784  $ abs1( beta( j )-bi( j, j ) ) /
785  $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
786  $ / ulp
787  IF( j.LT.mplusn ) THEN
788  IF( ai( j+1, j ).NE.zero ) THEN
789  ilabad = .true.
790  result( 5 ) = ulpinv
791  END IF
792  END IF
793  IF( j.GT.1 ) THEN
794  IF( ai( j, j-1 ).NE.zero ) THEN
795  ilabad = .true.
796  result( 5 ) = ulpinv
797  END IF
798  END IF
799  temp1 = max( temp1, temp2 )
800  IF( ilabad ) THEN
801  WRITE( nout, fmt = 9997 )j, mplusn, nptknt
802  END IF
803  110 CONTINUE
804  result( 6 ) = temp1
805 *
806 * Test (7) (if sorting worked) <--------- need to be checked.
807 *
808  ntest = 7
809  result( 7 ) = zero
810  IF( linfo.EQ.mplusn+3 )
811  $ result( 7 ) = ulpinv
812 *
813 * Test (8): compare the estimated value of DIF and its true value.
814 *
815  ntest = 8
816  result( 8 ) = zero
817  IF( difest( 2 ).EQ.zero ) THEN
818  IF( diftru.GT.abnrm*ulp )
819  $ result( 8 ) = ulpinv
820  ELSE IF( diftru.EQ.zero ) THEN
821  IF( difest( 2 ).GT.abnrm*ulp )
822  $ result( 8 ) = ulpinv
823  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
824  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
825  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
826  END IF
827 *
828 * Test (9)
829 *
830  ntest = 9
831  result( 9 ) = zero
832  IF( linfo.EQ.( mplusn+2 ) ) THEN
833  IF( diftru.GT.abnrm*ulp )
834  $ result( 9 ) = ulpinv
835  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
836  $ result( 9 ) = ulpinv
837  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
838  $ result( 9 ) = ulpinv
839  END IF
840 *
841 * Test (10): compare the estimated value of PL and it true value.
842 *
843  ntest = 10
844  result( 10 ) = zero
845  IF( pl( 1 ).EQ.zero ) THEN
846  IF( pltru.GT.abnrm*ulp )
847  $ result( 10 ) = ulpinv
848  ELSE IF( pltru.EQ.zero ) THEN
849  IF( pl( 1 ).GT.abnrm*ulp )
850  $ result( 10 ) = ulpinv
851  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
852  $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
853  result( 10 ) = ulpinv
854  END IF
855 *
856  ntestt = ntestt + ntest
857 *
858 * Print out tests which fail.
859 *
860  DO 120 j = 1, ntest
861  IF( result( j ).GE.thresh ) THEN
862 *
863 * If this is the first test to fail,
864 * print a header to the data file.
865 *
866  IF( nerrs.EQ.0 ) THEN
867  WRITE( nout, fmt = 9996 )'CGX'
868 *
869 * Matrix types
870 *
871  WRITE( nout, fmt = 9995 )
872 *
873 * Tests performed
874 *
875  WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
876  $ ( '''', i = 1, 4 )
877 *
878  END IF
879  nerrs = nerrs + 1
880  IF( result( j ).LT.10000.0 ) THEN
881  WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
882  ELSE
883  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
884  END IF
885  END IF
886 *
887  120 CONTINUE
888 *
889  130 CONTINUE
890  go to 80
891  140 CONTINUE
892 *
893  150 CONTINUE
894 *
895 * Summary
896 *
897  CALL alasvm( 'CGX', nout, nerrs, ntestt, 0 )
898 *
899  work( 1 ) = maxwrk
900 *
901  RETURN
902 *
903  9999 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
904  $ i6, ', JTYPE=', i6, ')' )
905 *
906  9998 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
907  $ i6, ', Input Example #', i2, ')' )
908 *
909  9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', i6, '.',
910  $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
911 *
912  9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
913  $ ' problem driver' )
914 *
915  9995 FORMAT( 'Input Example' )
916 *
917  9994 FORMAT( ' Matrix types: ', /
918  $ ' 1: A is a block diagonal matrix of Jordan blocks ',
919  $ 'and B is the identity ', / ' matrix, ',
920  $ / ' 2: A and B are upper triangular matrices, ',
921  $ / ' 3: A and B are as type 2, but each second diagonal ',
922  $ 'block in A_11 and ', /
923  $ ' each third diaongal block in A_22 are 2x2 blocks,',
924  $ / ' 4: A and B are block diagonal matrices, ',
925  $ / ' 5: (A,B) has potentially close or common ',
926  $ 'eigenvalues.', / )
927 *
928  9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
929  $ 'Q and Z are ', a, ',', / 19x,
930  $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
931  $ / ' 1 = | A - Q S Z', a,
932  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
933  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
934  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
935  $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
936  $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
937  $ ' and diagonals of (S,T)', /
938  $ ' 7 = 1/ULP if SDIM is not the correct number of ',
939  $ 'selected eigenvalues', /
940  $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
941  $ 'DIFTRU/DIFEST > 10*THRESH',
942  $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
943  $ 'when reordering fails', /
944  $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
945  $ 'PLTRU/PLEST > THRESH', /
946  $ ' ( Test 10 is only for input examples )', / )
947  9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
948  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
949  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
950  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
951  9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
952  $ ' result ', i2, ' is', 0p, f8.2 )
953  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
954  $ ' result ', i2, ' is', 1p, e10.3 )
955 *
956 * End of CDRGSX
957 *
958  END