LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
cget24.f
Go to the documentation of this file.
1 *> \brief \b CGET24
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 CGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
12 * H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
13 * RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
14 * LWORK, RWORK, BWORK, INFO )
15 *
16 * .. Scalar Arguments ..
17 * LOGICAL COMP
18 * INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
19 * $ NSLCT
20 * REAL RCDEIN, RCDVIN, THRESH
21 * ..
22 * .. Array Arguments ..
23 * LOGICAL BWORK( * )
24 * INTEGER ISEED( 4 ), ISLCT( * )
25 * REAL RESULT( 17 ), RWORK( * )
26 * COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
27 * $ VS( LDVS, * ), VS1( LDVS, * ), W( * ),
28 * $ WORK( * ), WT( * ), WTMP( * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> CGET24 checks the nonsymmetric eigenvalue (Schur form) problem
38 *> expert driver CGEESX.
39 *>
40 *> If COMP = .FALSE., the first 13 of the following tests will be
41 *> be performed on the input matrix A, and also tests 14 and 15
42 *> if LWORK is sufficiently large.
43 *> If COMP = .TRUE., all 17 test will be performed.
44 *>
45 *> (1) 0 if T is in Schur form, 1/ulp otherwise
46 *> (no sorting of eigenvalues)
47 *>
48 *> (2) | A - VS T VS' | / ( n |A| ulp )
49 *>
50 *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
51 *> form (no sorting of eigenvalues).
52 *>
53 *> (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
54 *>
55 *> (4) 0 if W are eigenvalues of T
56 *> 1/ulp otherwise
57 *> (no sorting of eigenvalues)
58 *>
59 *> (5) 0 if T(with VS) = T(without VS),
60 *> 1/ulp otherwise
61 *> (no sorting of eigenvalues)
62 *>
63 *> (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
64 *> 1/ulp otherwise
65 *> (no sorting of eigenvalues)
66 *>
67 *> (7) 0 if T is in Schur form, 1/ulp otherwise
68 *> (with sorting of eigenvalues)
69 *>
70 *> (8) | A - VS T VS' | / ( n |A| ulp )
71 *>
72 *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
73 *> form (with sorting of eigenvalues).
74 *>
75 *> (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
76 *>
77 *> (10) 0 if W are eigenvalues of T
78 *> 1/ulp otherwise
79 *> If workspace sufficient, also compare W with and
80 *> without reciprocal condition numbers
81 *> (with sorting of eigenvalues)
82 *>
83 *> (11) 0 if T(with VS) = T(without VS),
84 *> 1/ulp otherwise
85 *> If workspace sufficient, also compare T with and without
86 *> reciprocal condition numbers
87 *> (with sorting of eigenvalues)
88 *>
89 *> (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
90 *> 1/ulp otherwise
91 *> If workspace sufficient, also compare VS with and without
92 *> reciprocal condition numbers
93 *> (with sorting of eigenvalues)
94 *>
95 *> (13) if sorting worked and SDIM is the number of
96 *> eigenvalues which were SELECTed
97 *> If workspace sufficient, also compare SDIM with and
98 *> without reciprocal condition numbers
99 *>
100 *> (14) if RCONDE the same no matter if VS and/or RCONDV computed
101 *>
102 *> (15) if RCONDV the same no matter if VS and/or RCONDE computed
103 *>
104 *> (16) |RCONDE - RCDEIN| / cond(RCONDE)
105 *>
106 *> RCONDE is the reciprocal average eigenvalue condition number
107 *> computed by CGEESX and RCDEIN (the precomputed true value)
108 *> is supplied as input. cond(RCONDE) is the condition number
109 *> of RCONDE, and takes errors in computing RCONDE into account,
110 *> so that the resulting quantity should be O(ULP). cond(RCONDE)
111 *> is essentially given by norm(A)/RCONDV.
112 *>
113 *> (17) |RCONDV - RCDVIN| / cond(RCONDV)
114 *>
115 *> RCONDV is the reciprocal right invariant subspace condition
116 *> number computed by CGEESX and RCDVIN (the precomputed true
117 *> value) is supplied as input. cond(RCONDV) is the condition
118 *> number of RCONDV, and takes errors in computing RCONDV into
119 *> account, so that the resulting quantity should be O(ULP).
120 *> cond(RCONDV) is essentially given by norm(A)/RCONDE.
121 *> \endverbatim
122 *
123 * Arguments:
124 * ==========
125 *
126 *> \param[in] COMP
127 *> \verbatim
128 *> COMP is LOGICAL
129 *> COMP describes which input tests to perform:
130 *> = .FALSE. if the computed condition numbers are not to
131 *> be tested against RCDVIN and RCDEIN
132 *> = .TRUE. if they are to be compared
133 *> \endverbatim
134 *>
135 *> \param[in] JTYPE
136 *> \verbatim
137 *> JTYPE is INTEGER
138 *> Type of input matrix. Used to label output if error occurs.
139 *> \endverbatim
140 *>
141 *> \param[in] ISEED
142 *> \verbatim
143 *> ISEED is INTEGER array, dimension (4)
144 *> If COMP = .FALSE., the random number generator seed
145 *> used to produce matrix.
146 *> If COMP = .TRUE., ISEED(1) = the number of the example.
147 *> Used to label output if error occurs.
148 *> \endverbatim
149 *>
150 *> \param[in] THRESH
151 *> \verbatim
152 *> THRESH is REAL
153 *> A test will count as "failed" if the "error", computed as
154 *> described above, exceeds THRESH. Note that the error
155 *> is scaled to be O(1), so THRESH should be a reasonably
156 *> small multiple of 1, e.g., 10 or 100. In particular,
157 *> it should not depend on the precision (single vs. double)
158 *> or the size of the matrix. It must be at least zero.
159 *> \endverbatim
160 *>
161 *> \param[in] NOUNIT
162 *> \verbatim
163 *> NOUNIT is INTEGER
164 *> The FORTRAN unit number for printing out error messages
165 *> (e.g., if a routine returns INFO not equal to 0.)
166 *> \endverbatim
167 *>
168 *> \param[in] N
169 *> \verbatim
170 *> N is INTEGER
171 *> The dimension of A. N must be at least 0.
172 *> \endverbatim
173 *>
174 *> \param[in,out] A
175 *> \verbatim
176 *> A is COMPLEX array, dimension (LDA, N)
177 *> Used to hold the matrix whose eigenvalues are to be
178 *> computed.
179 *> \endverbatim
180 *>
181 *> \param[in] LDA
182 *> \verbatim
183 *> LDA is INTEGER
184 *> The leading dimension of A, and H. LDA must be at
185 *> least 1 and at least N.
186 *> \endverbatim
187 *>
188 *> \param[out] H
189 *> \verbatim
190 *> H is COMPLEX array, dimension (LDA, N)
191 *> Another copy of the test matrix A, modified by CGEESX.
192 *> \endverbatim
193 *>
194 *> \param[out] HT
195 *> \verbatim
196 *> HT is COMPLEX array, dimension (LDA, N)
197 *> Yet another copy of the test matrix A, modified by CGEESX.
198 *> \endverbatim
199 *>
200 *> \param[out] W
201 *> \verbatim
202 *> W is COMPLEX array, dimension (N)
203 *> The computed eigenvalues of A.
204 *> \endverbatim
205 *>
206 *> \param[out] WT
207 *> \verbatim
208 *> WT is COMPLEX array, dimension (N)
209 *> Like W, this array contains the eigenvalues of A,
210 *> but those computed when CGEESX only computes a partial
211 *> eigendecomposition, i.e. not Schur vectors
212 *> \endverbatim
213 *>
214 *> \param[out] WTMP
215 *> \verbatim
216 *> WTMP is COMPLEX array, dimension (N)
217 *> Like W, this array contains the eigenvalues of A,
218 *> but sorted by increasing real or imaginary part.
219 *> \endverbatim
220 *>
221 *> \param[out] VS
222 *> \verbatim
223 *> VS is COMPLEX array, dimension (LDVS, N)
224 *> VS holds the computed Schur vectors.
225 *> \endverbatim
226 *>
227 *> \param[in] LDVS
228 *> \verbatim
229 *> LDVS is INTEGER
230 *> Leading dimension of VS. Must be at least max(1, N).
231 *> \endverbatim
232 *>
233 *> \param[out] VS1
234 *> \verbatim
235 *> VS1 is COMPLEX array, dimension (LDVS, N)
236 *> VS1 holds another copy of the computed Schur vectors.
237 *> \endverbatim
238 *>
239 *> \param[in] RCDEIN
240 *> \verbatim
241 *> RCDEIN is REAL
242 *> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
243 *> condition number for the average of selected eigenvalues.
244 *> \endverbatim
245 *>
246 *> \param[in] RCDVIN
247 *> \verbatim
248 *> RCDVIN is REAL
249 *> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
250 *> condition number for the selected right invariant subspace.
251 *> \endverbatim
252 *>
253 *> \param[in] NSLCT
254 *> \verbatim
255 *> NSLCT is INTEGER
256 *> When COMP = .TRUE. the number of selected eigenvalues
257 *> corresponding to the precomputed values RCDEIN and RCDVIN.
258 *> \endverbatim
259 *>
260 *> \param[in] ISLCT
261 *> \verbatim
262 *> ISLCT is INTEGER array, dimension (NSLCT)
263 *> When COMP = .TRUE. ISLCT selects the eigenvalues of the
264 *> input matrix corresponding to the precomputed values RCDEIN
265 *> and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
266 *> eigenvalue with the J-th largest real or imaginary part is
267 *> selected. The real part is used if ISRT = 0, and the
268 *> imaginary part if ISRT = 1.
269 *> Not referenced if COMP = .FALSE.
270 *> \endverbatim
271 *>
272 *> \param[in] ISRT
273 *> \verbatim
274 *> ISRT is INTEGER
275 *> When COMP = .TRUE., ISRT describes how ISLCT is used to
276 *> choose a subset of the spectrum.
277 *> Not referenced if COMP = .FALSE.
278 *> \endverbatim
279 *>
280 *> \param[out] RESULT
281 *> \verbatim
282 *> RESULT is REAL array, dimension (17)
283 *> The values computed by the 17 tests described above.
284 *> The values are currently limited to 1/ulp, to avoid
285 *> overflow.
286 *> \endverbatim
287 *>
288 *> \param[out] WORK
289 *> \verbatim
290 *> WORK is COMPLEX array, dimension (2*N*N)
291 *> \endverbatim
292 *>
293 *> \param[in] LWORK
294 *> \verbatim
295 *> LWORK is INTEGER
296 *> The number of entries in WORK to be passed to CGEESX. This
297 *> must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to
298 *> be performed.
299 *> \endverbatim
300 *>
301 *> \param[out] RWORK
302 *> \verbatim
303 *> RWORK is REAL array, dimension (N)
304 *> \endverbatim
305 *>
306 *> \param[out] BWORK
307 *> \verbatim
308 *> BWORK is LOGICAL array, dimension (N)
309 *> \endverbatim
310 *>
311 *> \param[out] INFO
312 *> \verbatim
313 *> INFO is INTEGER
314 *> If 0, successful exit.
315 *> If <0, input parameter -INFO had an incorrect value.
316 *> If >0, CGEESX returned an error code, the absolute
317 *> value of which is returned.
318 *> \endverbatim
319 *
320 * Authors:
321 * ========
322 *
323 *> \author Univ. of Tennessee
324 *> \author Univ. of California Berkeley
325 *> \author Univ. of Colorado Denver
326 *> \author NAG Ltd.
327 *
328 *> \ingroup complex_eig
329 *
330 * =====================================================================
331  SUBROUTINE cget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
332  $ H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
333  $ RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
334  $ LWORK, RWORK, BWORK, INFO )
335 *
336 * -- LAPACK test routine --
337 * -- LAPACK is a software package provided by Univ. of Tennessee, --
338 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
339 *
340 * .. Scalar Arguments ..
341  LOGICAL COMP
342  INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
343  $ NSLCT
344  REAL RCDEIN, RCDVIN, THRESH
345 * ..
346 * .. Array Arguments ..
347  LOGICAL BWORK( * )
348  INTEGER ISEED( 4 ), ISLCT( * )
349  REAL RESULT( 17 ), RWORK( * )
350  COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
351  $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
352  $ work( * ), wt( * ), wtmp( * )
353 * ..
354 *
355 * =====================================================================
356 *
357 * .. Parameters ..
358  COMPLEX CZERO, CONE
359  PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
360  $ cone = ( 1.0e+0, 0.0e+0 ) )
361  REAL ZERO, ONE
362  parameter( zero = 0.0e+0, one = 1.0e+0 )
363  REAL EPSIN
364  parameter( epsin = 5.9605e-8 )
365 * ..
366 * .. Local Scalars ..
367  CHARACTER SORT
368  INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
369  $ SDIM, SDIM1
370  REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
371  $ smlnum, tol, tolin, ulp, ulpinv, v, vricmp,
372  $ vrimin, wnorm
373  COMPLEX CTMP
374 * ..
375 * .. Local Arrays ..
376  INTEGER IPNT( 20 )
377 * ..
378 * .. External Functions ..
379  LOGICAL CSLECT
380  REAL CLANGE, SLAMCH
381  EXTERNAL CSLECT, CLANGE, SLAMCH
382 * ..
383 * .. External Subroutines ..
384  EXTERNAL ccopy, cgeesx, cgemm, clacpy, cunt01, xerbla
385 * ..
386 * .. Intrinsic Functions ..
387  INTRINSIC abs, aimag, max, min, real
388 * ..
389 * .. Arrays in Common ..
390  LOGICAL SELVAL( 20 )
391  REAL SELWI( 20 ), SELWR( 20 )
392 * ..
393 * .. Scalars in Common ..
394  INTEGER SELDIM, SELOPT
395 * ..
396 * .. Common blocks ..
397  COMMON / sslct / selopt, seldim, selval, selwr, selwi
398 * ..
399 * .. Executable Statements ..
400 *
401 * Check for errors
402 *
403  info = 0
404  IF( thresh.LT.zero ) THEN
405  info = -3
406  ELSE IF( nounit.LE.0 ) THEN
407  info = -5
408  ELSE IF( n.LT.0 ) THEN
409  info = -6
410  ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
411  info = -8
412  ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
413  info = -15
414  ELSE IF( lwork.LT.2*n ) THEN
415  info = -24
416  END IF
417 *
418  IF( info.NE.0 ) THEN
419  CALL xerbla( 'CGET24', -info )
420  RETURN
421  END IF
422 *
423 * Quick return if nothing to do
424 *
425  DO 10 i = 1, 17
426  result( i ) = -one
427  10 CONTINUE
428 *
429  IF( n.EQ.0 )
430  $ RETURN
431 *
432 * Important constants
433 *
434  smlnum = slamch( 'Safe minimum' )
435  ulp = slamch( 'Precision' )
436  ulpinv = one / ulp
437 *
438 * Perform tests (1)-(13)
439 *
440  selopt = 0
441  DO 90 isort = 0, 1
442  IF( isort.EQ.0 ) THEN
443  sort = 'N'
444  rsub = 0
445  ELSE
446  sort = 'S'
447  rsub = 6
448  END IF
449 *
450 * Compute Schur form and Schur vectors, and test them
451 *
452  CALL clacpy( 'F', n, n, a, lda, h, lda )
453  CALL cgeesx( 'V', sort, cslect, 'N', n, h, lda, sdim, w, vs,
454  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
455  $ iinfo )
456  IF( iinfo.NE.0 ) THEN
457  result( 1+rsub ) = ulpinv
458  IF( jtype.NE.22 ) THEN
459  WRITE( nounit, fmt = 9998 )'CGEESX1', iinfo, n, jtype,
460  $ iseed
461  ELSE
462  WRITE( nounit, fmt = 9999 )'CGEESX1', iinfo, n,
463  $ iseed( 1 )
464  END IF
465  info = abs( iinfo )
466  RETURN
467  END IF
468  IF( isort.EQ.0 ) THEN
469  CALL ccopy( n, w, 1, wtmp, 1 )
470  END IF
471 *
472 * Do Test (1) or Test (7)
473 *
474  result( 1+rsub ) = zero
475  DO 30 j = 1, n - 1
476  DO 20 i = j + 1, n
477  IF( h( i, j ).NE.czero )
478  $ result( 1+rsub ) = ulpinv
479  20 CONTINUE
480  30 CONTINUE
481 *
482 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
483 *
484 * Copy A to VS1, used as workspace
485 *
486  CALL clacpy( ' ', n, n, a, lda, vs1, ldvs )
487 *
488 * Compute Q*H and store in HT.
489 *
490  CALL cgemm( 'No transpose', 'No transpose', n, n, n, cone, vs,
491  $ ldvs, h, lda, czero, ht, lda )
492 *
493 * Compute A - Q*H*Q'
494 *
495  CALL cgemm( 'No transpose', 'Conjugate transpose', n, n, n,
496  $ -cone, ht, lda, vs, ldvs, cone, vs1, ldvs )
497 *
498  anorm = max( clange( '1', n, n, a, lda, rwork ), smlnum )
499  wnorm = clange( '1', n, n, vs1, ldvs, rwork )
500 *
501  IF( anorm.GT.wnorm ) THEN
502  result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
503  ELSE
504  IF( anorm.LT.one ) THEN
505  result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
506  $ ( n*ulp )
507  ELSE
508  result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
509  $ ( n*ulp )
510  END IF
511  END IF
512 *
513 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
514 *
515  CALL cunt01( 'Columns', n, n, vs, ldvs, work, lwork, rwork,
516  $ result( 3+rsub ) )
517 *
518 * Do Test (4) or Test (10)
519 *
520  result( 4+rsub ) = zero
521  DO 40 i = 1, n
522  IF( h( i, i ).NE.w( i ) )
523  $ result( 4+rsub ) = ulpinv
524  40 CONTINUE
525 *
526 * Do Test (5) or Test (11)
527 *
528  CALL clacpy( 'F', n, n, a, lda, ht, lda )
529  CALL cgeesx( 'N', sort, cslect, 'N', n, ht, lda, sdim, wt, vs,
530  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
531  $ iinfo )
532  IF( iinfo.NE.0 ) THEN
533  result( 5+rsub ) = ulpinv
534  IF( jtype.NE.22 ) THEN
535  WRITE( nounit, fmt = 9998 )'CGEESX2', iinfo, n, jtype,
536  $ iseed
537  ELSE
538  WRITE( nounit, fmt = 9999 )'CGEESX2', iinfo, n,
539  $ iseed( 1 )
540  END IF
541  info = abs( iinfo )
542  GO TO 220
543  END IF
544 *
545  result( 5+rsub ) = zero
546  DO 60 j = 1, n
547  DO 50 i = 1, n
548  IF( h( i, j ).NE.ht( i, j ) )
549  $ result( 5+rsub ) = ulpinv
550  50 CONTINUE
551  60 CONTINUE
552 *
553 * Do Test (6) or Test (12)
554 *
555  result( 6+rsub ) = zero
556  DO 70 i = 1, n
557  IF( w( i ).NE.wt( i ) )
558  $ result( 6+rsub ) = ulpinv
559  70 CONTINUE
560 *
561 * Do Test (13)
562 *
563  IF( isort.EQ.1 ) THEN
564  result( 13 ) = zero
565  knteig = 0
566  DO 80 i = 1, n
567  IF( cslect( w( i ) ) )
568  $ knteig = knteig + 1
569  IF( i.LT.n ) THEN
570  IF( cslect( w( i+1 ) ) .AND.
571  $ ( .NOT.cslect( w( i ) ) ) )result( 13 ) = ulpinv
572  END IF
573  80 CONTINUE
574  IF( sdim.NE.knteig )
575  $ result( 13 ) = ulpinv
576  END IF
577 *
578  90 CONTINUE
579 *
580 * If there is enough workspace, perform tests (14) and (15)
581 * as well as (10) through (13)
582 *
583  IF( lwork.GE.( n*( n+1 ) ) / 2 ) THEN
584 *
585 * Compute both RCONDE and RCONDV with VS
586 *
587  sort = 'S'
588  result( 14 ) = zero
589  result( 15 ) = zero
590  CALL clacpy( 'F', n, n, a, lda, ht, lda )
591  CALL cgeesx( 'V', sort, cslect, 'B', n, ht, lda, sdim1, wt,
592  $ vs1, ldvs, rconde, rcondv, work, lwork, rwork,
593  $ bwork, iinfo )
594  IF( iinfo.NE.0 ) THEN
595  result( 14 ) = ulpinv
596  result( 15 ) = ulpinv
597  IF( jtype.NE.22 ) THEN
598  WRITE( nounit, fmt = 9998 )'CGEESX3', iinfo, n, jtype,
599  $ iseed
600  ELSE
601  WRITE( nounit, fmt = 9999 )'CGEESX3', iinfo, n,
602  $ iseed( 1 )
603  END IF
604  info = abs( iinfo )
605  GO TO 220
606  END IF
607 *
608 * Perform tests (10), (11), (12), and (13)
609 *
610  DO 110 i = 1, n
611  IF( w( i ).NE.wt( i ) )
612  $ result( 10 ) = ulpinv
613  DO 100 j = 1, n
614  IF( h( i, j ).NE.ht( i, j ) )
615  $ result( 11 ) = ulpinv
616  IF( vs( i, j ).NE.vs1( i, j ) )
617  $ result( 12 ) = ulpinv
618  100 CONTINUE
619  110 CONTINUE
620  IF( sdim.NE.sdim1 )
621  $ result( 13 ) = ulpinv
622 *
623 * Compute both RCONDE and RCONDV without VS, and compare
624 *
625  CALL clacpy( 'F', n, n, a, lda, ht, lda )
626  CALL cgeesx( 'N', sort, cslect, 'B', n, ht, lda, sdim1, wt,
627  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
628  $ bwork, iinfo )
629  IF( iinfo.NE.0 ) THEN
630  result( 14 ) = ulpinv
631  result( 15 ) = ulpinv
632  IF( jtype.NE.22 ) THEN
633  WRITE( nounit, fmt = 9998 )'CGEESX4', iinfo, n, jtype,
634  $ iseed
635  ELSE
636  WRITE( nounit, fmt = 9999 )'CGEESX4', iinfo, n,
637  $ iseed( 1 )
638  END IF
639  info = abs( iinfo )
640  GO TO 220
641  END IF
642 *
643 * Perform tests (14) and (15)
644 *
645  IF( rcnde1.NE.rconde )
646  $ result( 14 ) = ulpinv
647  IF( rcndv1.NE.rcondv )
648  $ result( 15 ) = ulpinv
649 *
650 * Perform tests (10), (11), (12), and (13)
651 *
652  DO 130 i = 1, n
653  IF( w( i ).NE.wt( i ) )
654  $ result( 10 ) = ulpinv
655  DO 120 j = 1, n
656  IF( h( i, j ).NE.ht( i, j ) )
657  $ result( 11 ) = ulpinv
658  IF( vs( i, j ).NE.vs1( i, j ) )
659  $ result( 12 ) = ulpinv
660  120 CONTINUE
661  130 CONTINUE
662  IF( sdim.NE.sdim1 )
663  $ result( 13 ) = ulpinv
664 *
665 * Compute RCONDE with VS, and compare
666 *
667  CALL clacpy( 'F', n, n, a, lda, ht, lda )
668  CALL cgeesx( 'V', sort, cslect, 'E', n, ht, lda, sdim1, wt,
669  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
670  $ bwork, iinfo )
671  IF( iinfo.NE.0 ) THEN
672  result( 14 ) = ulpinv
673  IF( jtype.NE.22 ) THEN
674  WRITE( nounit, fmt = 9998 )'CGEESX5', iinfo, n, jtype,
675  $ iseed
676  ELSE
677  WRITE( nounit, fmt = 9999 )'CGEESX5', iinfo, n,
678  $ iseed( 1 )
679  END IF
680  info = abs( iinfo )
681  GO TO 220
682  END IF
683 *
684 * Perform test (14)
685 *
686  IF( rcnde1.NE.rconde )
687  $ result( 14 ) = ulpinv
688 *
689 * Perform tests (10), (11), (12), and (13)
690 *
691  DO 150 i = 1, n
692  IF( w( i ).NE.wt( i ) )
693  $ result( 10 ) = ulpinv
694  DO 140 j = 1, n
695  IF( h( i, j ).NE.ht( i, j ) )
696  $ result( 11 ) = ulpinv
697  IF( vs( i, j ).NE.vs1( i, j ) )
698  $ result( 12 ) = ulpinv
699  140 CONTINUE
700  150 CONTINUE
701  IF( sdim.NE.sdim1 )
702  $ result( 13 ) = ulpinv
703 *
704 * Compute RCONDE without VS, and compare
705 *
706  CALL clacpy( 'F', n, n, a, lda, ht, lda )
707  CALL cgeesx( 'N', sort, cslect, 'E', n, ht, lda, sdim1, wt,
708  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
709  $ bwork, iinfo )
710  IF( iinfo.NE.0 ) THEN
711  result( 14 ) = ulpinv
712  IF( jtype.NE.22 ) THEN
713  WRITE( nounit, fmt = 9998 )'CGEESX6', iinfo, n, jtype,
714  $ iseed
715  ELSE
716  WRITE( nounit, fmt = 9999 )'CGEESX6', iinfo, n,
717  $ iseed( 1 )
718  END IF
719  info = abs( iinfo )
720  GO TO 220
721  END IF
722 *
723 * Perform test (14)
724 *
725  IF( rcnde1.NE.rconde )
726  $ result( 14 ) = ulpinv
727 *
728 * Perform tests (10), (11), (12), and (13)
729 *
730  DO 170 i = 1, n
731  IF( w( i ).NE.wt( i ) )
732  $ result( 10 ) = ulpinv
733  DO 160 j = 1, n
734  IF( h( i, j ).NE.ht( i, j ) )
735  $ result( 11 ) = ulpinv
736  IF( vs( i, j ).NE.vs1( i, j ) )
737  $ result( 12 ) = ulpinv
738  160 CONTINUE
739  170 CONTINUE
740  IF( sdim.NE.sdim1 )
741  $ result( 13 ) = ulpinv
742 *
743 * Compute RCONDV with VS, and compare
744 *
745  CALL clacpy( 'F', n, n, a, lda, ht, lda )
746  CALL cgeesx( 'V', sort, cslect, 'V', n, ht, lda, sdim1, wt,
747  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
748  $ bwork, iinfo )
749  IF( iinfo.NE.0 ) THEN
750  result( 15 ) = ulpinv
751  IF( jtype.NE.22 ) THEN
752  WRITE( nounit, fmt = 9998 )'CGEESX7', iinfo, n, jtype,
753  $ iseed
754  ELSE
755  WRITE( nounit, fmt = 9999 )'CGEESX7', iinfo, n,
756  $ iseed( 1 )
757  END IF
758  info = abs( iinfo )
759  GO TO 220
760  END IF
761 *
762 * Perform test (15)
763 *
764  IF( rcndv1.NE.rcondv )
765  $ result( 15 ) = ulpinv
766 *
767 * Perform tests (10), (11), (12), and (13)
768 *
769  DO 190 i = 1, n
770  IF( w( i ).NE.wt( i ) )
771  $ result( 10 ) = ulpinv
772  DO 180 j = 1, n
773  IF( h( i, j ).NE.ht( i, j ) )
774  $ result( 11 ) = ulpinv
775  IF( vs( i, j ).NE.vs1( i, j ) )
776  $ result( 12 ) = ulpinv
777  180 CONTINUE
778  190 CONTINUE
779  IF( sdim.NE.sdim1 )
780  $ result( 13 ) = ulpinv
781 *
782 * Compute RCONDV without VS, and compare
783 *
784  CALL clacpy( 'F', n, n, a, lda, ht, lda )
785  CALL cgeesx( 'N', sort, cslect, 'V', n, ht, lda, sdim1, wt,
786  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
787  $ bwork, iinfo )
788  IF( iinfo.NE.0 ) THEN
789  result( 15 ) = ulpinv
790  IF( jtype.NE.22 ) THEN
791  WRITE( nounit, fmt = 9998 )'CGEESX8', iinfo, n, jtype,
792  $ iseed
793  ELSE
794  WRITE( nounit, fmt = 9999 )'CGEESX8', iinfo, n,
795  $ iseed( 1 )
796  END IF
797  info = abs( iinfo )
798  GO TO 220
799  END IF
800 *
801 * Perform test (15)
802 *
803  IF( rcndv1.NE.rcondv )
804  $ result( 15 ) = ulpinv
805 *
806 * Perform tests (10), (11), (12), and (13)
807 *
808  DO 210 i = 1, n
809  IF( w( i ).NE.wt( i ) )
810  $ result( 10 ) = ulpinv
811  DO 200 j = 1, n
812  IF( h( i, j ).NE.ht( i, j ) )
813  $ result( 11 ) = ulpinv
814  IF( vs( i, j ).NE.vs1( i, j ) )
815  $ result( 12 ) = ulpinv
816  200 CONTINUE
817  210 CONTINUE
818  IF( sdim.NE.sdim1 )
819  $ result( 13 ) = ulpinv
820 *
821  END IF
822 *
823  220 CONTINUE
824 *
825 * If there are precomputed reciprocal condition numbers, compare
826 * computed values with them.
827 *
828  IF( comp ) THEN
829 *
830 * First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
831 * the logical function CSLECT selects the eigenvalues specified
832 * by NSLCT, ISLCT and ISRT.
833 *
834  seldim = n
835  selopt = 1
836  eps = max( ulp, epsin )
837  DO 230 i = 1, n
838  ipnt( i ) = i
839  selval( i ) = .false.
840  selwr( i ) = real( wtmp( i ) )
841  selwi( i ) = aimag( wtmp( i ) )
842  230 CONTINUE
843  DO 250 i = 1, n - 1
844  kmin = i
845  IF( isrt.EQ.0 ) THEN
846  vrimin = real( wtmp( i ) )
847  ELSE
848  vrimin = aimag( wtmp( i ) )
849  END IF
850  DO 240 j = i + 1, n
851  IF( isrt.EQ.0 ) THEN
852  vricmp = real( wtmp( j ) )
853  ELSE
854  vricmp = aimag( wtmp( j ) )
855  END IF
856  IF( vricmp.LT.vrimin ) THEN
857  kmin = j
858  vrimin = vricmp
859  END IF
860  240 CONTINUE
861  ctmp = wtmp( kmin )
862  wtmp( kmin ) = wtmp( i )
863  wtmp( i ) = ctmp
864  itmp = ipnt( i )
865  ipnt( i ) = ipnt( kmin )
866  ipnt( kmin ) = itmp
867  250 CONTINUE
868  DO 260 i = 1, nslct
869  selval( ipnt( islct( i ) ) ) = .true.
870  260 CONTINUE
871 *
872 * Compute condition numbers
873 *
874  CALL clacpy( 'F', n, n, a, lda, ht, lda )
875  CALL cgeesx( 'N', 'S', cslect, 'B', n, ht, lda, sdim1, wt, vs1,
876  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
877  $ iinfo )
878  IF( iinfo.NE.0 ) THEN
879  result( 16 ) = ulpinv
880  result( 17 ) = ulpinv
881  WRITE( nounit, fmt = 9999 )'CGEESX9', iinfo, n, iseed( 1 )
882  info = abs( iinfo )
883  GO TO 270
884  END IF
885 *
886 * Compare condition number for average of selected eigenvalues
887 * taking its condition number into account
888 *
889  anorm = clange( '1', n, n, a, lda, rwork )
890  v = max( real( n )*eps*anorm, smlnum )
891  IF( anorm.EQ.zero )
892  $ v = one
893  IF( v.GT.rcondv ) THEN
894  tol = one
895  ELSE
896  tol = v / rcondv
897  END IF
898  IF( v.GT.rcdvin ) THEN
899  tolin = one
900  ELSE
901  tolin = v / rcdvin
902  END IF
903  tol = max( tol, smlnum / eps )
904  tolin = max( tolin, smlnum / eps )
905  IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
906  result( 16 ) = ulpinv
907  ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
908  result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
909  ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
910  result( 16 ) = ulpinv
911  ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
912  result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
913  ELSE
914  result( 16 ) = one
915  END IF
916 *
917 * Compare condition numbers for right invariant subspace
918 * taking its condition number into account
919 *
920  IF( v.GT.rcondv*rconde ) THEN
921  tol = rcondv
922  ELSE
923  tol = v / rconde
924  END IF
925  IF( v.GT.rcdvin*rcdein ) THEN
926  tolin = rcdvin
927  ELSE
928  tolin = v / rcdein
929  END IF
930  tol = max( tol, smlnum / eps )
931  tolin = max( tolin, smlnum / eps )
932  IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
933  result( 17 ) = ulpinv
934  ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
935  result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
936  ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
937  result( 17 ) = ulpinv
938  ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
939  result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
940  ELSE
941  result( 17 ) = one
942  END IF
943 *
944  270 CONTINUE
945 *
946  END IF
947 *
948  9999 FORMAT( ' CGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
949  $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
950  9998 FORMAT( ' CGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
951  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
952 *
953  RETURN
954 *
955 * End of CGET24
956 *
957  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine cunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
CUNT01
Definition: cunt01.f:126
subroutine cget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK, LWORK, RWORK, BWORK, INFO)
CGET24
Definition: cget24.f:335
subroutine cgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, INFO)
CGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: cgeesx.f:239
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