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