LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine cggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: cggesx.f:332
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvd.f:216
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
Definition: clakf2.f:107
subroutine cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
Definition: cget51.f:156
subroutine clatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
CLATM5
Definition: clatm5.f:270
subroutine cdrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
CDRGSX
Definition: cdrgsx.f:351