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