LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
zdrgsx.f
Go to the documentation of this file.
1 *> \brief \b ZDRGSX
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 ZDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
12 * BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK,
13 * RWORK, 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 RWORK( * ), S( * )
24 * COMPLEX*16 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 *> ZDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
37 *> problem expert driver ZGGESX.
38 *>
39 *> ZGGES 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 ZDRGSX is called with NSIZE > 0, five (5) types of built-in
55 *> matrix pairs are used to test the routine ZGGESX.
56 *>
57 *> When ZDRGSX is called with NSIZE = 0, it reads in test matrix data
58 *> to test ZGGESX.
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 ZGGESX, 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 ZGESVD) 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 DLATM5).
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 DGGESX.
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 ZLAKF2
189 *> \endverbatim
190 *>
191 *> \param[in] THRESH
192 *> \verbatim
193 *> THRESH is DOUBLE PRECISION
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*16 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*16 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*16 array, dimension (LDA, NSIZE)
241 *> Copy of A, modified by ZGGESX.
242 *> \endverbatim
243 *>
244 *> \param[out] BI
245 *> \verbatim
246 *> BI is COMPLEX*16 array, dimension (LDA, NSIZE)
247 *> Copy of B, modified by ZGGESX.
248 *> \endverbatim
249 *>
250 *> \param[out] Z
251 *> \verbatim
252 *> Z is COMPLEX*16 array, dimension (LDA, NSIZE)
253 *> Z holds the left Schur vectors computed by ZGGESX.
254 *> \endverbatim
255 *>
256 *> \param[out] Q
257 *> \verbatim
258 *> Q is COMPLEX*16 array, dimension (LDA, NSIZE)
259 *> Q holds the right Schur vectors computed by ZGGESX.
260 *> \endverbatim
261 *>
262 *> \param[out] ALPHA
263 *> \verbatim
264 *> ALPHA is COMPLEX*16 array, dimension (NSIZE)
265 *> \endverbatim
266 *>
267 *> \param[out] BETA
268 *> \verbatim
269 *> BETA is COMPLEX*16 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*16 array, dimension (LDC, LDC)
277 *> Store the matrix generated by subroutine ZLAKF2, 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 DOUBLE PRECISION array, dimension (LDC)
291 *> Singular values of C
292 *> \endverbatim
293 *>
294 *> \param[out] WORK
295 *> \verbatim
296 *> WORK is COMPLEX*16 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 DOUBLE PRECISION 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 complex16_eig
346 *
347 * =====================================================================
348  SUBROUTINE zdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
349  $ bi, z, q, alpha, beta, c, ldc, s, work, lwork,
350  $ 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  DOUBLE PRECISION THRESH
361 * ..
362 * .. Array Arguments ..
363  LOGICAL BWORK( * )
364  INTEGER IWORK( * )
365  DOUBLE PRECISION RWORK( * ), S( * )
366  COMPLEX*16 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  DOUBLE PRECISION ZERO, ONE, TEN
376  parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
377  COMPLEX*16 CZERO
378  parameter( czero = ( 0.0d+0, 0.0d+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  DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
387  $ temp2, thrsh2, ulp, ulpinv, weight
388  COMPLEX*16 X
389 * ..
390 * .. Local Arrays ..
391  DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
392 * ..
393 * .. External Functions ..
394  LOGICAL ZLCTSX
395  INTEGER ILAENV
396  DOUBLE PRECISION DLAMCH, ZLANGE
397  EXTERNAL zlctsx, ilaenv, dlamch, zlange
398 * ..
399 * .. External Subroutines ..
400  EXTERNAL alasvm, dlabad, xerbla, zgesvd, zget51, zggesx,
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, dble, dimag, max, sqrt
412 * ..
413 * .. Statement Functions ..
414  DOUBLE PRECISION ABS1
415 * ..
416 * .. Statement Function definitions ..
417  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
418 * ..
419 * .. Executable Statements ..
420 *
421 * Check for errors
422 *
423  info = 0
424  IF( nsize.LT.0 ) THEN
425  info = -1
426  ELSE IF( thresh.LT.zero ) THEN
427  info = -2
428  ELSE IF( nin.LE.0 ) THEN
429  info = -3
430  ELSE IF( nout.LE.0 ) THEN
431  info = -4
432  ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
433  info = -6
434  ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
435  info = -15
436  ELSE IF( liwork.LT.nsize+2 ) THEN
437  info = -21
438  END IF
439 *
440 * Compute workspace
441 * (Note: Comments in the code beginning "Workspace:" describe the
442 * minimal amount of workspace needed at that point in the code,
443 * as well as the preferred amount for good performance.
444 * NB refers to the optimal block size for the immediately
445 * following subroutine, as returned by ILAENV.)
446 *
447  minwrk = 1
448  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
449  minwrk = 3*nsize*nsize / 2
450 *
451 * workspace for cggesx
452 *
453  maxwrk = nsize*( 1+ilaenv( 1, 'ZGEQRF', ' ', nsize, 1, nsize,
454  $ 0 ) )
455  maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1, 'ZUNGQR', ' ',
456  $ nsize, 1, nsize, -1 ) ) )
457 *
458 * workspace for zgesvd
459 *
460  bdspac = 3*nsize*nsize / 2
461  maxwrk = max( maxwrk, nsize*nsize*
462  $ ( 1+ilaenv( 1, 'ZGEBRD', ' ', 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 = -18
473 *
474  IF( info.NE.0 ) THEN
475  CALL xerbla( 'ZDRGSX', -info )
476  RETURN
477  END IF
478 *
479 * Important constants
480 *
481  ulp = dlamch( 'P' )
482  ulpinv = one / ulp
483  smlnum = dlamch( 'S' ) / ulp
484  bignum = one / smlnum
485  CALL dlabad( 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 ZGGESX, 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 zlaset( 'Full', mplusn, mplusn, czero, czero, ai,
519  $ lda )
520  CALL zlaset( 'Full', mplusn, mplusn, czero, czero, bi,
521  $ lda )
522 *
523  CALL zlatm5( 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 ZLCTSX
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 zlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
544  CALL zlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
545 *
546  CALL zggesx( 'V', 'V', 'S', zlctsx, sense, mplusn, ai,
547  $ lda, bi, lda, mm, alpha, beta, q, lda, z,
548  $ lda, pl, difest, work, lwork, rwork,
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 )'ZGGESX', linfo, mplusn,
554  $ prtype
555  info = linfo
556  GO TO 30
557  END IF
558 *
559 * Compute the norm(A, B)
560 *
561  CALL zlacpy( 'Full', mplusn, mplusn, ai, lda, work,
562  $ mplusn )
563  CALL zlacpy( 'Full', mplusn, mplusn, bi, lda,
564  $ work( mplusn*mplusn+1 ), mplusn )
565  abnrm = zlange( 'Fro', mplusn, 2*mplusn, work, mplusn,
566  $ rwork )
567 *
568 * Do tests (1) to (4)
569 *
570  result( 2 ) = zero
571  CALL zget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
572  $ lda, work, rwork, result( 1 ) )
573  CALL zget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
574  $ lda, work, rwork, result( 2 ) )
575  CALL zget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
576  $ lda, work, rwork, result( 3 ) )
577  CALL zget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
578  $ lda, work, rwork, result( 4 ) )
579  ntest = 4
580 *
581 * Do tests (5) and (6): check Schur form of A and
582 * compare eigenvalues with diagonals.
583 *
584  temp1 = zero
585  result( 5 ) = zero
586  result( 6 ) = zero
587 *
588  DO 10 j = 1, mplusn
589  ilabad = .false.
590  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
591  $ max( smlnum, abs1( alpha( j ) ),
592  $ abs1( ai( j, j ) ) )+
593  $ abs1( beta( j )-bi( j, j ) ) /
594  $ max( smlnum, abs1( beta( j ) ),
595  $ abs1( 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  temp1 = max( temp1, temp2 )
609  IF( ilabad ) THEN
610  WRITE( nout, fmt = 9997 )j, mplusn, prtype
611  END IF
612  10 CONTINUE
613  result( 6 ) = temp1
614  ntest = ntest + 2
615 *
616 * Test (7) (if sorting worked)
617 *
618  result( 7 ) = zero
619  IF( linfo.EQ.mplusn+3 ) THEN
620  result( 7 ) = ulpinv
621  ELSE IF( mm.NE.n ) THEN
622  result( 7 ) = ulpinv
623  END IF
624  ntest = ntest + 1
625 *
626 * Test (8): compare the estimated value DIF and its
627 * value. first, compute the exact DIF.
628 *
629  result( 8 ) = zero
630  mn2 = mm*( mplusn-mm )*2
631  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
632 *
633 * Note: for either following two cases, there are
634 * almost same number of test cases fail the test.
635 *
636  CALL zlakf2( mm, mplusn-mm, ai, lda,
637  $ ai( mm+1, mm+1 ), bi,
638  $ bi( mm+1, mm+1 ), c, ldc )
639 *
640  CALL zgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
641  $ 1, work( 2 ), 1, work( 3 ), lwork-2,
642  $ rwork, info )
643  diftru = s( mn2 )
644 *
645  IF( difest( 2 ).EQ.zero ) THEN
646  IF( diftru.GT.abnrm*ulp )
647  $ result( 8 ) = ulpinv
648  ELSE IF( diftru.EQ.zero ) THEN
649  IF( difest( 2 ).GT.abnrm*ulp )
650  $ result( 8 ) = ulpinv
651  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
652  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
653  result( 8 ) = max( diftru / difest( 2 ),
654  $ difest( 2 ) / diftru )
655  END IF
656  ntest = ntest + 1
657  END IF
658 *
659 * Test (9)
660 *
661  result( 9 ) = zero
662  IF( linfo.EQ.( mplusn+2 ) ) THEN
663  IF( diftru.GT.abnrm*ulp )
664  $ result( 9 ) = ulpinv
665  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
666  $ result( 9 ) = ulpinv
667  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
668  $ result( 9 ) = ulpinv
669  ntest = ntest + 1
670  END IF
671 *
672  ntestt = ntestt + ntest
673 *
674 * Print out tests which fail.
675 *
676  DO 20 j = 1, 9
677  IF( result( j ).GE.thresh ) THEN
678 *
679 * If this is the first test to fail,
680 * print a header to the data file.
681 *
682  IF( nerrs.EQ.0 ) THEN
683  WRITE( nout, fmt = 9996 )'ZGX'
684 *
685 * Matrix types
686 *
687  WRITE( nout, fmt = 9994 )
688 *
689 * Tests performed
690 *
691  WRITE( nout, fmt = 9993 )'unitary', '''',
692  $ 'transpose', ( '''', i = 1, 4 )
693 *
694  END IF
695  nerrs = nerrs + 1
696  IF( result( j ).LT.10000.0d0 ) THEN
697  WRITE( nout, fmt = 9992 )mplusn, prtype,
698  $ weight, m, j, result( j )
699  ELSE
700  WRITE( nout, fmt = 9991 )mplusn, prtype,
701  $ weight, m, j, result( j )
702  END IF
703  END IF
704  20 CONTINUE
705 *
706  30 CONTINUE
707  40 CONTINUE
708  50 CONTINUE
709  60 CONTINUE
710 *
711  GO TO 150
712 *
713  70 CONTINUE
714 *
715 * Read in data from file to check accuracy of condition estimation
716 * Read input data until N=0
717 *
718  nptknt = 0
719 *
720  80 CONTINUE
721  READ( nin, fmt = *, end = 140 )mplusn
722  IF( mplusn.EQ.0 )
723  $ GO TO 140
724  READ( nin, fmt = *, end = 140 )n
725  DO 90 i = 1, mplusn
726  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
727  90 CONTINUE
728  DO 100 i = 1, mplusn
729  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
730  100 CONTINUE
731  READ( nin, fmt = * )pltru, diftru
732 *
733  nptknt = nptknt + 1
734  fs = .true.
735  k = 0
736  m = mplusn - n
737 *
738  CALL zlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
739  CALL zlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
740 *
741 * Compute the Schur factorization while swaping the
742 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
743 *
744  CALL zggesx( 'V', 'V', 'S', zlctsx, 'B', mplusn, ai, lda, bi, lda,
745  $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
746  $ lwork, rwork, iwork, liwork, bwork, linfo )
747 *
748  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
749  result( 1 ) = ulpinv
750  WRITE( nout, fmt = 9998 )'ZGGESX', linfo, mplusn, nptknt
751  GO TO 130
752  END IF
753 *
754 * Compute the norm(A, B)
755 * (should this be norm of (A,B) or (AI,BI)?)
756 *
757  CALL zlacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
758  CALL zlacpy( 'Full', mplusn, mplusn, bi, lda,
759  $ work( mplusn*mplusn+1 ), mplusn )
760  abnrm = zlange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
761 *
762 * Do tests (1) to (4)
763 *
764  CALL zget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
765  $ rwork, result( 1 ) )
766  CALL zget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
767  $ rwork, result( 2 ) )
768  CALL zget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
769  $ rwork, result( 3 ) )
770  CALL zget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
771  $ rwork, result( 4 ) )
772 *
773 * Do tests (5) and (6): check Schur form of A and compare
774 * eigenvalues with diagonals.
775 *
776  ntest = 6
777  temp1 = zero
778  result( 5 ) = zero
779  result( 6 ) = zero
780 *
781  DO 110 j = 1, mplusn
782  ilabad = .false.
783  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
784  $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
785  $ abs1( beta( j )-bi( j, j ) ) /
786  $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
787  $ / ulp
788  IF( j.LT.mplusn ) THEN
789  IF( ai( j+1, j ).NE.zero ) THEN
790  ilabad = .true.
791  result( 5 ) = ulpinv
792  END IF
793  END IF
794  IF( j.GT.1 ) THEN
795  IF( ai( j, j-1 ).NE.zero ) THEN
796  ilabad = .true.
797  result( 5 ) = ulpinv
798  END IF
799  END IF
800  temp1 = max( temp1, temp2 )
801  IF( ilabad ) THEN
802  WRITE( nout, fmt = 9997 )j, mplusn, nptknt
803  END IF
804  110 CONTINUE
805  result( 6 ) = temp1
806 *
807 * Test (7) (if sorting worked) <--------- need to be checked.
808 *
809  ntest = 7
810  result( 7 ) = zero
811  IF( linfo.EQ.mplusn+3 )
812  $ result( 7 ) = ulpinv
813 *
814 * Test (8): compare the estimated value of DIF and its true value.
815 *
816  ntest = 8
817  result( 8 ) = zero
818  IF( difest( 2 ).EQ.zero ) THEN
819  IF( diftru.GT.abnrm*ulp )
820  $ result( 8 ) = ulpinv
821  ELSE IF( diftru.EQ.zero ) THEN
822  IF( difest( 2 ).GT.abnrm*ulp )
823  $ result( 8 ) = ulpinv
824  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
825  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
826  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
827  END IF
828 *
829 * Test (9)
830 *
831  ntest = 9
832  result( 9 ) = zero
833  IF( linfo.EQ.( mplusn+2 ) ) THEN
834  IF( diftru.GT.abnrm*ulp )
835  $ result( 9 ) = ulpinv
836  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
837  $ result( 9 ) = ulpinv
838  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
839  $ result( 9 ) = ulpinv
840  END IF
841 *
842 * Test (10): compare the estimated value of PL and it true value.
843 *
844  ntest = 10
845  result( 10 ) = zero
846  IF( pl( 1 ).EQ.zero ) THEN
847  IF( pltru.GT.abnrm*ulp )
848  $ result( 10 ) = ulpinv
849  ELSE IF( pltru.EQ.zero ) THEN
850  IF( pl( 1 ).GT.abnrm*ulp )
851  $ result( 10 ) = ulpinv
852  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
853  $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
854  result( 10 ) = ulpinv
855  END IF
856 *
857  ntestt = ntestt + ntest
858 *
859 * Print out tests which fail.
860 *
861  DO 120 j = 1, ntest
862  IF( result( j ).GE.thresh ) THEN
863 *
864 * If this is the first test to fail,
865 * print a header to the data file.
866 *
867  IF( nerrs.EQ.0 ) THEN
868  WRITE( nout, fmt = 9996 )'ZGX'
869 *
870 * Matrix types
871 *
872  WRITE( nout, fmt = 9995 )
873 *
874 * Tests performed
875 *
876  WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
877  $ ( '''', i = 1, 4 )
878 *
879  END IF
880  nerrs = nerrs + 1
881  IF( result( j ).LT.10000.0d0 ) THEN
882  WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
883  ELSE
884  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
885  END IF
886  END IF
887 *
888  120 CONTINUE
889 *
890  130 CONTINUE
891  GO TO 80
892  140 CONTINUE
893 *
894  150 CONTINUE
895 *
896 * Summary
897 *
898  CALL alasvm( 'ZGX', nout, nerrs, ntestt, 0 )
899 *
900  work( 1 ) = maxwrk
901 *
902  RETURN
903 *
904  9999 FORMAT( ' ZDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
905  $ i6, ', JTYPE=', i6, ')' )
906 *
907  9998 FORMAT( ' ZDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
908  $ i6, ', Input Example #', i2, ')' )
909 *
910  9997 FORMAT( ' ZDRGSX: S not in Schur form at eigenvalue ', i6, '.',
911  $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
912 *
913  9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
914  $ ' problem driver' )
915 *
916  9995 FORMAT( 'Input Example' )
917 *
918  9994 FORMAT( ' Matrix types: ', /
919  $ ' 1: A is a block diagonal matrix of Jordan blocks ',
920  $ 'and B is the identity ', / ' matrix, ',
921  $ / ' 2: A and B are upper triangular matrices, ',
922  $ / ' 3: A and B are as type 2, but each second diagonal ',
923  $ 'block in A_11 and ', /
924  $ ' each third diaongal block in A_22 are 2x2 blocks,',
925  $ / ' 4: A and B are block diagonal matrices, ',
926  $ / ' 5: (A,B) has potentially close or common ',
927  $ 'eigenvalues.', / )
928 *
929  9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
930  $ 'Q and Z are ', a, ',', / 19x,
931  $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
932  $ / ' 1 = | A - Q S Z', a,
933  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
934  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
935  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
936  $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
937  $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
938  $ ' and diagonals of (S,T)', /
939  $ ' 7 = 1/ULP if SDIM is not the correct number of ',
940  $ 'selected eigenvalues', /
941  $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
942  $ 'DIFTRU/DIFEST > 10*THRESH',
943  $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
944  $ 'when reordering fails', /
945  $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
946  $ 'PLTRU/PLEST > THRESH', /
947  $ ' ( Test 10 is only for input examples )', / )
948  9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
949  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
950  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
951  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, d10.3 )
952  9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
953  $ ' result ', i2, ' is', 0p, f8.2 )
954  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
955  $ ' result ', i2, ' is', 1p, d10.3 )
956 *
957 * End of ZDRGSX
958 *
959  END
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zggesx(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)
ZGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: zggesx.f:332
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine zget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
ZGET51
Definition: zget51.f:156
subroutine zlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
ZLAKF2
Definition: zlakf2.f:107
subroutine zgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: zgesvd.f:216
subroutine zdrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
ZDRGSX
Definition: zdrgsx.f:351
subroutine zlatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
ZLATM5
Definition: zlatm5.f:270