LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sdrgsx.f
Go to the documentation of this file.
1 *> \brief \b SDRGSX
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 SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
12 * AI, 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 * REAL THRESH
19 * ..
20 * .. Array Arguments ..
21 * LOGICAL BWORK( * )
22 * INTEGER IWORK( * )
23 * REAL 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 *> SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
36 *> problem expert driver SGGESX.
37 *>
38 *> SGGESX 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 SDRGSX is called with NSIZE > 0, five (5) types of built-in
57 *> matrix pairs are used to test the routine SGGESX.
58 *>
59 *> When SDRGSX is called with NSIZE = 0, it reads in test matrix data
60 *> to test SGGESX.
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 SGGESX, 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 SGESVD) 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 SLATM5).
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 SGGESX.
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 SLAKF2
199 *> \endverbatim
200 *>
201 *> \param[in] THRESH
202 *> \verbatim
203 *> THRESH is REAL
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 REAL 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 REAL 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 REAL array, dimension (LDA, NSIZE)
251 *> Copy of A, modified by SGGESX.
252 *> \endverbatim
253 *>
254 *> \param[out] BI
255 *> \verbatim
256 *> BI is REAL array, dimension (LDA, NSIZE)
257 *> Copy of B, modified by SGGESX.
258 *> \endverbatim
259 *>
260 *> \param[out] Z
261 *> \verbatim
262 *> Z is REAL array, dimension (LDA, NSIZE)
263 *> Z holds the left Schur vectors computed by SGGESX.
264 *> \endverbatim
265 *>
266 *> \param[out] Q
267 *> \verbatim
268 *> Q is REAL array, dimension (LDA, NSIZE)
269 *> Q holds the right Schur vectors computed by SGGESX.
270 *> \endverbatim
271 *>
272 *> \param[out] ALPHAR
273 *> \verbatim
274 *> ALPHAR is REAL array, dimension (NSIZE)
275 *> \endverbatim
276 *>
277 *> \param[out] ALPHAI
278 *> \verbatim
279 *> ALPHAI is REAL array, dimension (NSIZE)
280 *> \endverbatim
281 *>
282 *> \param[out] BETA
283 *> \verbatim
284 *> BETA is REAL array, dimension (NSIZE)
285 *> \verbatim
286 *> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
287 *> \endverbatim
288 *>
289 *> \param[out] C
290 *> \verbatim
291 *> C is REAL array, dimension (LDC, LDC)
292 *> Store the matrix generated by subroutine SLAKF2, 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 REAL array, dimension (LDC)
306 *> Singular values of C
307 *> \endverbatim
308 *>
309 *> \param[out] WORK
310 *> \verbatim
311 *> WORK is REAL 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 single_eig
356 *
357 * =====================================================================
358  SUBROUTINE sdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
359  $ ai, 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  REAL thresh
371 * ..
372 * .. Array Arguments ..
373  LOGICAL bwork( * )
374  INTEGER iwork( * )
375  REAL 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  REAL zero, one, ten
385  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+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  REAL abnrm, bignum, diftru, pltru, smlnum, temp1,
394  $ temp2, thrsh2, ulp, ulpinv, weight
395 * ..
396 * .. Local Arrays ..
397  REAL difest( 2 ), pl( 2 ), result( 10 )
398 * ..
399 * .. External Functions ..
400  LOGICAL slctsx
401  INTEGER ilaenv
402  REAL slamch, slange
403  EXTERNAL slctsx, ilaenv, slamch, slange
404 * ..
405 * .. External Subroutines ..
406  EXTERNAL alasvm, sgesvd, sget51, sget53, sggesx, slabad,
407  $ slacpy, slakf2, slaset, slatm5, xerbla
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 c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
449  minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
450 *
451 * workspace for sggesx
452 *
453  maxwrk = 9*( nsize+1 ) + nsize*
454  $ ilaenv( 1, 'SGEQRF', ' ', nsize, 1, nsize, 0 )
455  maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
456  $ ilaenv( 1, 'SORGQR', ' ', nsize, 1, nsize, -1 ) )
457 *
458 * workspace for sgesvd
459 *
460  bdspac = 5*nsize*nsize / 2
461  maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
462  $ ilaenv( 1, 'SGEBRD', ' ', 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 = -19
473 *
474  IF( info.NE.0 ) THEN
475  CALL xerbla( 'SDRGSX', -info )
476  return
477  END IF
478 *
479 * Important constants
480 *
481  ulp = slamch( 'P' )
482  ulpinv = one / ulp
483  smlnum = slamch( 'S' ) / ulp
484  bignum = one / smlnum
485  CALL slabad( 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 SGGESX, 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 slaset( 'Full', mplusn, mplusn, zero, zero, ai,
519  $ lda )
520  CALL slaset( 'Full', mplusn, mplusn, zero, zero, bi,
521  $ lda )
522 *
523  CALL slatm5( 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 SLCTSX
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 slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
544  CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
545 *
546  CALL sggesx( 'V', 'V', 'S', slctsx, sense, mplusn, ai,
547  $ lda, bi, lda, mm, alphar, alphai, beta,
548  $ q, lda, z, lda, pl, difest, work, lwork,
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 )'SGGESX', linfo, mplusn,
554  $ prtype
555  info = linfo
556  go to 30
557  END IF
558 *
559 * Compute the norm(A, B)
560 *
561  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work,
562  $ mplusn )
563  CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
564  $ work( mplusn*mplusn+1 ), mplusn )
565  abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn,
566  $ work )
567 *
568 * Do tests (1) to (4)
569 *
570  CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
571  $ lda, work, result( 1 ) )
572  CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
573  $ lda, work, result( 2 ) )
574  CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
575  $ lda, work, result( 3 ) )
576  CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
577  $ lda, work, 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  IF( alphai( j ).EQ.zero ) THEN
590  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
591  $ max( smlnum, abs( alphar( j ) ),
592  $ abs( ai( j, j ) ) )+
593  $ abs( beta( j )-bi( j, j ) ) /
594  $ max( smlnum, abs( beta( j ) ),
595  $ abs( 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  ELSE
609  IF( alphai( j ).GT.zero ) THEN
610  i1 = j
611  ELSE
612  i1 = j - 1
613  END IF
614  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
615  ilabad = .true.
616  ELSE IF( i1.LT.mplusn-1 ) THEN
617  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
618  ilabad = .true.
619  result( 5 ) = ulpinv
620  END IF
621  ELSE IF( i1.GT.1 ) THEN
622  IF( ai( i1, i1-1 ).NE.zero ) THEN
623  ilabad = .true.
624  result( 5 ) = ulpinv
625  END IF
626  END IF
627  IF( .NOT.ilabad ) THEN
628  CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ),
629  $ lda, beta( j ), alphar( j ),
630  $ alphai( j ), temp2, iinfo )
631  IF( iinfo.GE.3 ) THEN
632  WRITE( nout, fmt = 9997 )iinfo, j,
633  $ mplusn, prtype
634  info = abs( iinfo )
635  END IF
636  ELSE
637  temp2 = ulpinv
638  END IF
639  END IF
640  temp1 = max( temp1, temp2 )
641  IF( ilabad ) THEN
642  WRITE( nout, fmt = 9996 )j, mplusn, prtype
643  END IF
644  10 continue
645  result( 6 ) = temp1
646  ntest = ntest + 2
647 *
648 * Test (7) (if sorting worked)
649 *
650  result( 7 ) = zero
651  IF( linfo.EQ.mplusn+3 ) THEN
652  result( 7 ) = ulpinv
653  ELSE IF( mm.NE.n ) THEN
654  result( 7 ) = ulpinv
655  END IF
656  ntest = ntest + 1
657 *
658 * Test (8): compare the estimated value DIF and its
659 * value. first, compute the exact DIF.
660 *
661  result( 8 ) = zero
662  mn2 = mm*( mplusn-mm )*2
663  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
664 *
665 * Note: for either following two causes, there are
666 * almost same number of test cases fail the test.
667 *
668  CALL slakf2( mm, mplusn-mm, ai, lda,
669  $ ai( mm+1, mm+1 ), bi,
670  $ bi( mm+1, mm+1 ), c, ldc )
671 *
672  CALL sgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
673  $ 1, work( 2 ), 1, work( 3 ), lwork-2,
674  $ info )
675  diftru = s( mn2 )
676 *
677  IF( difest( 2 ).EQ.zero ) THEN
678  IF( diftru.GT.abnrm*ulp )
679  $ result( 8 ) = ulpinv
680  ELSE IF( diftru.EQ.zero ) THEN
681  IF( difest( 2 ).GT.abnrm*ulp )
682  $ result( 8 ) = ulpinv
683  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
684  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
685  result( 8 ) = max( diftru / difest( 2 ),
686  $ difest( 2 ) / diftru )
687  END IF
688  ntest = ntest + 1
689  END IF
690 *
691 * Test (9)
692 *
693  result( 9 ) = zero
694  IF( linfo.EQ.( mplusn+2 ) ) THEN
695  IF( diftru.GT.abnrm*ulp )
696  $ result( 9 ) = ulpinv
697  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
698  $ result( 9 ) = ulpinv
699  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
700  $ result( 9 ) = ulpinv
701  ntest = ntest + 1
702  END IF
703 *
704  ntestt = ntestt + ntest
705 *
706 * Print out tests which fail.
707 *
708  DO 20 j = 1, 9
709  IF( result( j ).GE.thresh ) THEN
710 *
711 * If this is the first test to fail,
712 * print a header to the data file.
713 *
714  IF( nerrs.EQ.0 ) THEN
715  WRITE( nout, fmt = 9995 )'SGX'
716 *
717 * Matrix types
718 *
719  WRITE( nout, fmt = 9993 )
720 *
721 * Tests performed
722 *
723  WRITE( nout, fmt = 9992 )'orthogonal', '''',
724  $ 'transpose', ( '''', i = 1, 4 )
725 *
726  END IF
727  nerrs = nerrs + 1
728  IF( result( j ).LT.10000.0 ) THEN
729  WRITE( nout, fmt = 9991 )mplusn, prtype,
730  $ weight, m, j, result( j )
731  ELSE
732  WRITE( nout, fmt = 9990 )mplusn, prtype,
733  $ weight, m, j, result( j )
734  END IF
735  END IF
736  20 continue
737 *
738  30 continue
739  40 continue
740  50 continue
741  60 continue
742 *
743  go to 150
744 *
745  70 continue
746 *
747 * Read in data from file to check accuracy of condition estimation
748 * Read input data until N=0
749 *
750  nptknt = 0
751 *
752  80 continue
753  READ( nin, fmt = *, END = 140 )mplusn
754  IF( mplusn.EQ.0 )
755  $ go to 140
756  READ( nin, fmt = *, END = 140 )n
757  DO 90 i = 1, mplusn
758  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
759  90 continue
760  DO 100 i = 1, mplusn
761  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
762  100 continue
763  READ( nin, fmt = * )pltru, diftru
764 *
765  nptknt = nptknt + 1
766  fs = .true.
767  k = 0
768  m = mplusn - n
769 *
770  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
771  CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
772 *
773 * Compute the Schur factorization while swaping the
774 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
775 *
776  CALL sggesx( 'V', 'V', 'S', slctsx, 'B', mplusn, ai, lda, bi, lda,
777  $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
778  $ work, lwork, iwork, liwork, bwork, linfo )
779 *
780  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
781  result( 1 ) = ulpinv
782  WRITE( nout, fmt = 9998 )'SGGESX', linfo, mplusn, nptknt
783  go to 130
784  END IF
785 *
786 * Compute the norm(A, B)
787 * (should this be norm of (A,B) or (AI,BI)?)
788 *
789  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
790  CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
791  $ work( mplusn*mplusn+1 ), mplusn )
792  abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
793 *
794 * Do tests (1) to (4)
795 *
796  CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
797  $ result( 1 ) )
798  CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
799  $ result( 2 ) )
800  CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
801  $ result( 3 ) )
802  CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
803  $ result( 4 ) )
804 *
805 * Do tests (5) and (6): check Schur form of A and compare
806 * eigenvalues with diagonals.
807 *
808  ntest = 6
809  temp1 = zero
810  result( 5 ) = zero
811  result( 6 ) = zero
812 *
813  DO 110 j = 1, mplusn
814  ilabad = .false.
815  IF( alphai( j ).EQ.zero ) THEN
816  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
817  $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
818  $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
819  $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
820  $ / ulp
821  IF( j.LT.mplusn ) THEN
822  IF( ai( j+1, j ).NE.zero ) THEN
823  ilabad = .true.
824  result( 5 ) = ulpinv
825  END IF
826  END IF
827  IF( j.GT.1 ) THEN
828  IF( ai( j, j-1 ).NE.zero ) THEN
829  ilabad = .true.
830  result( 5 ) = ulpinv
831  END IF
832  END IF
833  ELSE
834  IF( alphai( j ).GT.zero ) THEN
835  i1 = j
836  ELSE
837  i1 = j - 1
838  END IF
839  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
840  ilabad = .true.
841  ELSE IF( i1.LT.mplusn-1 ) THEN
842  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
843  ilabad = .true.
844  result( 5 ) = ulpinv
845  END IF
846  ELSE IF( i1.GT.1 ) THEN
847  IF( ai( i1, i1-1 ).NE.zero ) THEN
848  ilabad = .true.
849  result( 5 ) = ulpinv
850  END IF
851  END IF
852  IF( .NOT.ilabad ) THEN
853  CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
854  $ beta( j ), alphar( j ), alphai( j ), temp2,
855  $ iinfo )
856  IF( iinfo.GE.3 ) THEN
857  WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
858  info = abs( iinfo )
859  END IF
860  ELSE
861  temp2 = ulpinv
862  END IF
863  END IF
864  temp1 = max( temp1, temp2 )
865  IF( ilabad ) THEN
866  WRITE( nout, fmt = 9996 )j, mplusn, nptknt
867  END IF
868  110 continue
869  result( 6 ) = temp1
870 *
871 * Test (7) (if sorting worked) <--------- need to be checked.
872 *
873  ntest = 7
874  result( 7 ) = zero
875  IF( linfo.EQ.mplusn+3 )
876  $ result( 7 ) = ulpinv
877 *
878 * Test (8): compare the estimated value of DIF and its true value.
879 *
880  ntest = 8
881  result( 8 ) = zero
882  IF( difest( 2 ).EQ.zero ) THEN
883  IF( diftru.GT.abnrm*ulp )
884  $ result( 8 ) = ulpinv
885  ELSE IF( diftru.EQ.zero ) THEN
886  IF( difest( 2 ).GT.abnrm*ulp )
887  $ result( 8 ) = ulpinv
888  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
889  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
890  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
891  END IF
892 *
893 * Test (9)
894 *
895  ntest = 9
896  result( 9 ) = zero
897  IF( linfo.EQ.( mplusn+2 ) ) THEN
898  IF( diftru.GT.abnrm*ulp )
899  $ result( 9 ) = ulpinv
900  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
901  $ result( 9 ) = ulpinv
902  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
903  $ result( 9 ) = ulpinv
904  END IF
905 *
906 * Test (10): compare the estimated value of PL and it true value.
907 *
908  ntest = 10
909  result( 10 ) = zero
910  IF( pl( 1 ).EQ.zero ) THEN
911  IF( pltru.GT.abnrm*ulp )
912  $ result( 10 ) = ulpinv
913  ELSE IF( pltru.EQ.zero ) THEN
914  IF( pl( 1 ).GT.abnrm*ulp )
915  $ result( 10 ) = ulpinv
916  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
917  $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
918  result( 10 ) = ulpinv
919  END IF
920 *
921  ntestt = ntestt + ntest
922 *
923 * Print out tests which fail.
924 *
925  DO 120 j = 1, ntest
926  IF( result( j ).GE.thresh ) THEN
927 *
928 * If this is the first test to fail,
929 * print a header to the data file.
930 *
931  IF( nerrs.EQ.0 ) THEN
932  WRITE( nout, fmt = 9995 )'SGX'
933 *
934 * Matrix types
935 *
936  WRITE( nout, fmt = 9994 )
937 *
938 * Tests performed
939 *
940  WRITE( nout, fmt = 9992 )'orthogonal', '''',
941  $ 'transpose', ( '''', i = 1, 4 )
942 *
943  END IF
944  nerrs = nerrs + 1
945  IF( result( j ).LT.10000.0 ) THEN
946  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
947  ELSE
948  WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
949  END IF
950  END IF
951 *
952  120 continue
953 *
954  130 continue
955  go to 80
956  140 continue
957 *
958  150 continue
959 *
960 * Summary
961 *
962  CALL alasvm( 'SGX', nout, nerrs, ntestt, 0 )
963 *
964  work( 1 ) = maxwrk
965 *
966  return
967 *
968  9999 format( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
969  $ i6, ', JTYPE=', i6, ')' )
970 *
971  9998 format( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
972  $ i6, ', Input Example #', i2, ')' )
973 *
974  9997 format( ' SDRGSX: SGET53 returned INFO=', i1, ' for eigenvalue ',
975  $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
976 *
977  9996 format( ' SDRGSX: S not in Schur form at eigenvalue ', i6, '.',
978  $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
979 *
980  9995 format( / 1x, a3, ' -- Real Expert Generalized Schur form',
981  $ ' problem driver' )
982 *
983  9994 format( 'Input Example' )
984 *
985  9993 format( ' Matrix types: ', /
986  $ ' 1: A is a block diagonal matrix of Jordan blocks ',
987  $ 'and B is the identity ', / ' matrix, ',
988  $ / ' 2: A and B are upper triangular matrices, ',
989  $ / ' 3: A and B are as type 2, but each second diagonal ',
990  $ 'block in A_11 and ', /
991  $ ' each third diaongal block in A_22 are 2x2 blocks,',
992  $ / ' 4: A and B are block diagonal matrices, ',
993  $ / ' 5: (A,B) has potentially close or common ',
994  $ 'eigenvalues.', / )
995 *
996  9992 format( / ' Tests performed: (S is Schur, T is triangular, ',
997  $ 'Q and Z are ', a, ',', / 19x,
998  $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
999  $ / ' 1 = | A - Q S Z', a,
1000  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1001  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
1002  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
1003  $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1004  $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1005  $ ' and diagonals of (S,T)', /
1006  $ ' 7 = 1/ULP if SDIM is not the correct number of ',
1007  $ 'selected eigenvalues', /
1008  $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1009  $ 'DIFTRU/DIFEST > 10*THRESH',
1010  $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1011  $ 'when reordering fails', /
1012  $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1013  $ 'PLTRU/PLEST > THRESH', /
1014  $ ' ( Test 10 is only for input examples )', / )
1015  9991 format( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1016  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1017  9990 format( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1018  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
1019  9989 format( ' Input example #', i2, ', matrix order=', i4, ',',
1020  $ ' result ', i2, ' is', 0p, f8.2 )
1021  9988 format( ' Input example #', i2, ', matrix order=', i4, ',',
1022  $ ' result ', i2, ' is', 1p, e10.3 )
1023 *
1024 * End of SDRGSX
1025 *
1026  END