LAPACK  3.6.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 *> \date November 2011
329 *
330 *> \ingroup complex_eig
331 *
332 * =====================================================================
333  SUBROUTINE cget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
334  $ h, ht, w, wt, wtmp, vs, ldvs, vs1, rcdein,
335  $ rcdvin, nslct, islct, isrt, result, work,
336  $ lwork, rwork, bwork, info )
337 *
338 * -- LAPACK test routine (version 3.4.0) --
339 * -- LAPACK is a software package provided by Univ. of Tennessee, --
340 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
341 * November 2011
342 *
343 * .. Scalar Arguments ..
344  LOGICAL COMP
345  INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
346  $ nslct
347  REAL RCDEIN, RCDVIN, THRESH
348 * ..
349 * .. Array Arguments ..
350  LOGICAL BWORK( * )
351  INTEGER ISEED( 4 ), ISLCT( * )
352  REAL RESULT( 17 ), RWORK( * )
353  COMPLEX A( lda, * ), H( lda, * ), HT( lda, * ),
354  $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
355  $ work( * ), wt( * ), wtmp( * )
356 * ..
357 *
358 * =====================================================================
359 *
360 * .. Parameters ..
361  COMPLEX CZERO, CONE
362  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
363  $ cone = ( 1.0e+0, 0.0e+0 ) )
364  REAL ZERO, ONE
365  parameter ( zero = 0.0e+0, one = 1.0e+0 )
366  REAL EPSIN
367  parameter ( epsin = 5.9605e-8 )
368 * ..
369 * .. Local Scalars ..
370  CHARACTER SORT
371  INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
372  $ sdim, sdim1
373  REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
374  $ smlnum, tol, tolin, ulp, ulpinv, v, vricmp,
375  $ vrimin, wnorm
376  COMPLEX CTMP
377 * ..
378 * .. Local Arrays ..
379  INTEGER IPNT( 20 )
380 * ..
381 * .. External Functions ..
382  LOGICAL CSLECT
383  REAL CLANGE, SLAMCH
384  EXTERNAL cslect, clange, slamch
385 * ..
386 * .. External Subroutines ..
387  EXTERNAL ccopy, cgeesx, cgemm, clacpy, cunt01, xerbla
388 * ..
389 * .. Intrinsic Functions ..
390  INTRINSIC abs, aimag, max, min, real
391 * ..
392 * .. Arrays in Common ..
393  LOGICAL SELVAL( 20 )
394  REAL SELWI( 20 ), SELWR( 20 )
395 * ..
396 * .. Scalars in Common ..
397  INTEGER SELDIM, SELOPT
398 * ..
399 * .. Common blocks ..
400  COMMON / sslct / selopt, seldim, selval, selwr, selwi
401 * ..
402 * .. Executable Statements ..
403 *
404 * Check for errors
405 *
406  info = 0
407  IF( thresh.LT.zero ) THEN
408  info = -3
409  ELSE IF( nounit.LE.0 ) THEN
410  info = -5
411  ELSE IF( n.LT.0 ) THEN
412  info = -6
413  ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
414  info = -8
415  ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
416  info = -15
417  ELSE IF( lwork.LT.2*n ) THEN
418  info = -24
419  END IF
420 *
421  IF( info.NE.0 ) THEN
422  CALL xerbla( 'CGET24', -info )
423  RETURN
424  END IF
425 *
426 * Quick return if nothing to do
427 *
428  DO 10 i = 1, 17
429  result( i ) = -one
430  10 CONTINUE
431 *
432  IF( n.EQ.0 )
433  $ RETURN
434 *
435 * Important constants
436 *
437  smlnum = slamch( 'Safe minimum' )
438  ulp = slamch( 'Precision' )
439  ulpinv = one / ulp
440 *
441 * Perform tests (1)-(13)
442 *
443  selopt = 0
444  DO 90 isort = 0, 1
445  IF( isort.EQ.0 ) THEN
446  sort = 'N'
447  rsub = 0
448  ELSE
449  sort = 'S'
450  rsub = 6
451  END IF
452 *
453 * Compute Schur form and Schur vectors, and test them
454 *
455  CALL clacpy( 'F', n, n, a, lda, h, lda )
456  CALL cgeesx( 'V', sort, cslect, 'N', n, h, lda, sdim, w, vs,
457  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
458  $ iinfo )
459  IF( iinfo.NE.0 ) THEN
460  result( 1+rsub ) = ulpinv
461  IF( jtype.NE.22 ) THEN
462  WRITE( nounit, fmt = 9998 )'CGEESX1', iinfo, n, jtype,
463  $ iseed
464  ELSE
465  WRITE( nounit, fmt = 9999 )'CGEESX1', iinfo, n,
466  $ iseed( 1 )
467  END IF
468  info = abs( iinfo )
469  RETURN
470  END IF
471  IF( isort.EQ.0 ) THEN
472  CALL ccopy( n, w, 1, wtmp, 1 )
473  END IF
474 *
475 * Do Test (1) or Test (7)
476 *
477  result( 1+rsub ) = zero
478  DO 30 j = 1, n - 1
479  DO 20 i = j + 1, n
480  IF( h( i, j ).NE.czero )
481  $ result( 1+rsub ) = ulpinv
482  20 CONTINUE
483  30 CONTINUE
484 *
485 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
486 *
487 * Copy A to VS1, used as workspace
488 *
489  CALL clacpy( ' ', n, n, a, lda, vs1, ldvs )
490 *
491 * Compute Q*H and store in HT.
492 *
493  CALL cgemm( 'No transpose', 'No transpose', n, n, n, cone, vs,
494  $ ldvs, h, lda, czero, ht, lda )
495 *
496 * Compute A - Q*H*Q'
497 *
498  CALL cgemm( 'No transpose', 'Conjugate transpose', n, n, n,
499  $ -cone, ht, lda, vs, ldvs, cone, vs1, ldvs )
500 *
501  anorm = max( clange( '1', n, n, a, lda, rwork ), smlnum )
502  wnorm = clange( '1', n, n, vs1, ldvs, rwork )
503 *
504  IF( anorm.GT.wnorm ) THEN
505  result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
506  ELSE
507  IF( anorm.LT.one ) THEN
508  result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
509  $ ( n*ulp )
510  ELSE
511  result( 2+rsub ) = min( wnorm / anorm, REAL( N ) ) /
512  $ ( n*ulp )
513  END IF
514  END IF
515 *
516 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
517 *
518  CALL cunt01( 'Columns', n, n, vs, ldvs, work, lwork, rwork,
519  $ result( 3+rsub ) )
520 *
521 * Do Test (4) or Test (10)
522 *
523  result( 4+rsub ) = zero
524  DO 40 i = 1, n
525  IF( h( i, i ).NE.w( i ) )
526  $ result( 4+rsub ) = ulpinv
527  40 CONTINUE
528 *
529 * Do Test (5) or Test (11)
530 *
531  CALL clacpy( 'F', n, n, a, lda, ht, lda )
532  CALL cgeesx( 'N', sort, cslect, 'N', n, ht, lda, sdim, wt, vs,
533  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
534  $ iinfo )
535  IF( iinfo.NE.0 ) THEN
536  result( 5+rsub ) = ulpinv
537  IF( jtype.NE.22 ) THEN
538  WRITE( nounit, fmt = 9998 )'CGEESX2', iinfo, n, jtype,
539  $ iseed
540  ELSE
541  WRITE( nounit, fmt = 9999 )'CGEESX2', iinfo, n,
542  $ iseed( 1 )
543  END IF
544  info = abs( iinfo )
545  GO TO 220
546  END IF
547 *
548  result( 5+rsub ) = zero
549  DO 60 j = 1, n
550  DO 50 i = 1, n
551  IF( h( i, j ).NE.ht( i, j ) )
552  $ result( 5+rsub ) = ulpinv
553  50 CONTINUE
554  60 CONTINUE
555 *
556 * Do Test (6) or Test (12)
557 *
558  result( 6+rsub ) = zero
559  DO 70 i = 1, n
560  IF( w( i ).NE.wt( i ) )
561  $ result( 6+rsub ) = ulpinv
562  70 CONTINUE
563 *
564 * Do Test (13)
565 *
566  IF( isort.EQ.1 ) THEN
567  result( 13 ) = zero
568  knteig = 0
569  DO 80 i = 1, n
570  IF( cslect( w( i ) ) )
571  $ knteig = knteig + 1
572  IF( i.LT.n ) THEN
573  IF( cslect( w( i+1 ) ) .AND.
574  $ ( .NOT.cslect( w( i ) ) ) )result( 13 ) = ulpinv
575  END IF
576  80 CONTINUE
577  IF( sdim.NE.knteig )
578  $ result( 13 ) = ulpinv
579  END IF
580 *
581  90 CONTINUE
582 *
583 * If there is enough workspace, perform tests (14) and (15)
584 * as well as (10) through (13)
585 *
586  IF( lwork.GE.( n*( n+1 ) ) / 2 ) THEN
587 *
588 * Compute both RCONDE and RCONDV with VS
589 *
590  sort = 'S'
591  result( 14 ) = zero
592  result( 15 ) = zero
593  CALL clacpy( 'F', n, n, a, lda, ht, lda )
594  CALL cgeesx( 'V', sort, cslect, 'B', n, ht, lda, sdim1, wt,
595  $ vs1, ldvs, rconde, rcondv, work, lwork, rwork,
596  $ bwork, iinfo )
597  IF( iinfo.NE.0 ) THEN
598  result( 14 ) = ulpinv
599  result( 15 ) = ulpinv
600  IF( jtype.NE.22 ) THEN
601  WRITE( nounit, fmt = 9998 )'CGEESX3', iinfo, n, jtype,
602  $ iseed
603  ELSE
604  WRITE( nounit, fmt = 9999 )'CGEESX3', iinfo, n,
605  $ iseed( 1 )
606  END IF
607  info = abs( iinfo )
608  GO TO 220
609  END IF
610 *
611 * Perform tests (10), (11), (12), and (13)
612 *
613  DO 110 i = 1, n
614  IF( w( i ).NE.wt( i ) )
615  $ result( 10 ) = ulpinv
616  DO 100 j = 1, n
617  IF( h( i, j ).NE.ht( i, j ) )
618  $ result( 11 ) = ulpinv
619  IF( vs( i, j ).NE.vs1( i, j ) )
620  $ result( 12 ) = ulpinv
621  100 CONTINUE
622  110 CONTINUE
623  IF( sdim.NE.sdim1 )
624  $ result( 13 ) = ulpinv
625 *
626 * Compute both RCONDE and RCONDV without VS, and compare
627 *
628  CALL clacpy( 'F', n, n, a, lda, ht, lda )
629  CALL cgeesx( 'N', sort, cslect, 'B', n, ht, lda, sdim1, wt,
630  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
631  $ bwork, iinfo )
632  IF( iinfo.NE.0 ) THEN
633  result( 14 ) = ulpinv
634  result( 15 ) = ulpinv
635  IF( jtype.NE.22 ) THEN
636  WRITE( nounit, fmt = 9998 )'CGEESX4', iinfo, n, jtype,
637  $ iseed
638  ELSE
639  WRITE( nounit, fmt = 9999 )'CGEESX4', iinfo, n,
640  $ iseed( 1 )
641  END IF
642  info = abs( iinfo )
643  GO TO 220
644  END IF
645 *
646 * Perform tests (14) and (15)
647 *
648  IF( rcnde1.NE.rconde )
649  $ result( 14 ) = ulpinv
650  IF( rcndv1.NE.rcondv )
651  $ result( 15 ) = ulpinv
652 *
653 * Perform tests (10), (11), (12), and (13)
654 *
655  DO 130 i = 1, n
656  IF( w( i ).NE.wt( i ) )
657  $ result( 10 ) = ulpinv
658  DO 120 j = 1, n
659  IF( h( i, j ).NE.ht( i, j ) )
660  $ result( 11 ) = ulpinv
661  IF( vs( i, j ).NE.vs1( i, j ) )
662  $ result( 12 ) = ulpinv
663  120 CONTINUE
664  130 CONTINUE
665  IF( sdim.NE.sdim1 )
666  $ result( 13 ) = ulpinv
667 *
668 * Compute RCONDE with VS, and compare
669 *
670  CALL clacpy( 'F', n, n, a, lda, ht, lda )
671  CALL cgeesx( 'V', sort, cslect, 'E', n, ht, lda, sdim1, wt,
672  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
673  $ bwork, iinfo )
674  IF( iinfo.NE.0 ) THEN
675  result( 14 ) = ulpinv
676  IF( jtype.NE.22 ) THEN
677  WRITE( nounit, fmt = 9998 )'CGEESX5', iinfo, n, jtype,
678  $ iseed
679  ELSE
680  WRITE( nounit, fmt = 9999 )'CGEESX5', iinfo, n,
681  $ iseed( 1 )
682  END IF
683  info = abs( iinfo )
684  GO TO 220
685  END IF
686 *
687 * Perform test (14)
688 *
689  IF( rcnde1.NE.rconde )
690  $ result( 14 ) = ulpinv
691 *
692 * Perform tests (10), (11), (12), and (13)
693 *
694  DO 150 i = 1, n
695  IF( w( i ).NE.wt( i ) )
696  $ result( 10 ) = ulpinv
697  DO 140 j = 1, n
698  IF( h( i, j ).NE.ht( i, j ) )
699  $ result( 11 ) = ulpinv
700  IF( vs( i, j ).NE.vs1( i, j ) )
701  $ result( 12 ) = ulpinv
702  140 CONTINUE
703  150 CONTINUE
704  IF( sdim.NE.sdim1 )
705  $ result( 13 ) = ulpinv
706 *
707 * Compute RCONDE without VS, and compare
708 *
709  CALL clacpy( 'F', n, n, a, lda, ht, lda )
710  CALL cgeesx( 'N', sort, cslect, 'E', n, ht, lda, sdim1, wt,
711  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
712  $ bwork, iinfo )
713  IF( iinfo.NE.0 ) THEN
714  result( 14 ) = ulpinv
715  IF( jtype.NE.22 ) THEN
716  WRITE( nounit, fmt = 9998 )'CGEESX6', iinfo, n, jtype,
717  $ iseed
718  ELSE
719  WRITE( nounit, fmt = 9999 )'CGEESX6', iinfo, n,
720  $ iseed( 1 )
721  END IF
722  info = abs( iinfo )
723  GO TO 220
724  END IF
725 *
726 * Perform test (14)
727 *
728  IF( rcnde1.NE.rconde )
729  $ result( 14 ) = ulpinv
730 *
731 * Perform tests (10), (11), (12), and (13)
732 *
733  DO 170 i = 1, n
734  IF( w( i ).NE.wt( i ) )
735  $ result( 10 ) = ulpinv
736  DO 160 j = 1, n
737  IF( h( i, j ).NE.ht( i, j ) )
738  $ result( 11 ) = ulpinv
739  IF( vs( i, j ).NE.vs1( i, j ) )
740  $ result( 12 ) = ulpinv
741  160 CONTINUE
742  170 CONTINUE
743  IF( sdim.NE.sdim1 )
744  $ result( 13 ) = ulpinv
745 *
746 * Compute RCONDV with VS, and compare
747 *
748  CALL clacpy( 'F', n, n, a, lda, ht, lda )
749  CALL cgeesx( 'V', sort, cslect, 'V', n, ht, lda, sdim1, wt,
750  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
751  $ bwork, iinfo )
752  IF( iinfo.NE.0 ) THEN
753  result( 15 ) = ulpinv
754  IF( jtype.NE.22 ) THEN
755  WRITE( nounit, fmt = 9998 )'CGEESX7', iinfo, n, jtype,
756  $ iseed
757  ELSE
758  WRITE( nounit, fmt = 9999 )'CGEESX7', iinfo, n,
759  $ iseed( 1 )
760  END IF
761  info = abs( iinfo )
762  GO TO 220
763  END IF
764 *
765 * Perform test (15)
766 *
767  IF( rcndv1.NE.rcondv )
768  $ result( 15 ) = ulpinv
769 *
770 * Perform tests (10), (11), (12), and (13)
771 *
772  DO 190 i = 1, n
773  IF( w( i ).NE.wt( i ) )
774  $ result( 10 ) = ulpinv
775  DO 180 j = 1, n
776  IF( h( i, j ).NE.ht( i, j ) )
777  $ result( 11 ) = ulpinv
778  IF( vs( i, j ).NE.vs1( i, j ) )
779  $ result( 12 ) = ulpinv
780  180 CONTINUE
781  190 CONTINUE
782  IF( sdim.NE.sdim1 )
783  $ result( 13 ) = ulpinv
784 *
785 * Compute RCONDV without VS, and compare
786 *
787  CALL clacpy( 'F', n, n, a, lda, ht, lda )
788  CALL cgeesx( 'N', sort, cslect, 'V', n, ht, lda, sdim1, wt,
789  $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
790  $ bwork, iinfo )
791  IF( iinfo.NE.0 ) THEN
792  result( 15 ) = ulpinv
793  IF( jtype.NE.22 ) THEN
794  WRITE( nounit, fmt = 9998 )'CGEESX8', iinfo, n, jtype,
795  $ iseed
796  ELSE
797  WRITE( nounit, fmt = 9999 )'CGEESX8', iinfo, n,
798  $ iseed( 1 )
799  END IF
800  info = abs( iinfo )
801  GO TO 220
802  END IF
803 *
804 * Perform test (15)
805 *
806  IF( rcndv1.NE.rcondv )
807  $ result( 15 ) = ulpinv
808 *
809 * Perform tests (10), (11), (12), and (13)
810 *
811  DO 210 i = 1, n
812  IF( w( i ).NE.wt( i ) )
813  $ result( 10 ) = ulpinv
814  DO 200 j = 1, n
815  IF( h( i, j ).NE.ht( i, j ) )
816  $ result( 11 ) = ulpinv
817  IF( vs( i, j ).NE.vs1( i, j ) )
818  $ result( 12 ) = ulpinv
819  200 CONTINUE
820  210 CONTINUE
821  IF( sdim.NE.sdim1 )
822  $ result( 13 ) = ulpinv
823 *
824  END IF
825 *
826  220 CONTINUE
827 *
828 * If there are precomputed reciprocal condition numbers, compare
829 * computed values with them.
830 *
831  IF( comp ) THEN
832 *
833 * First set up SELOPT, SELDIM, SELVAL, SELWR and SELWI so that
834 * the logical function CSLECT selects the eigenvalues specified
835 * by NSLCT, ISLCT and ISRT.
836 *
837  seldim = n
838  selopt = 1
839  eps = max( ulp, epsin )
840  DO 230 i = 1, n
841  ipnt( i ) = i
842  selval( i ) = .false.
843  selwr( i ) = REAL( WTMP( I ) )
844  selwi( i ) = aimag( wtmp( i ) )
845  230 CONTINUE
846  DO 250 i = 1, n - 1
847  kmin = i
848  IF( isrt.EQ.0 ) THEN
849  vrimin = REAL( WTMP( I ) )
850  ELSE
851  vrimin = aimag( wtmp( i ) )
852  END IF
853  DO 240 j = i + 1, n
854  IF( isrt.EQ.0 ) THEN
855  vricmp = REAL( WTMP( J ) )
856  ELSE
857  vricmp = aimag( wtmp( j ) )
858  END IF
859  IF( vricmp.LT.vrimin ) THEN
860  kmin = j
861  vrimin = vricmp
862  END IF
863  240 CONTINUE
864  ctmp = wtmp( kmin )
865  wtmp( kmin ) = wtmp( i )
866  wtmp( i ) = ctmp
867  itmp = ipnt( i )
868  ipnt( i ) = ipnt( kmin )
869  ipnt( kmin ) = itmp
870  250 CONTINUE
871  DO 260 i = 1, nslct
872  selval( ipnt( islct( i ) ) ) = .true.
873  260 CONTINUE
874 *
875 * Compute condition numbers
876 *
877  CALL clacpy( 'F', n, n, a, lda, ht, lda )
878  CALL cgeesx( 'N', 'S', cslect, 'B', n, ht, lda, sdim1, wt, vs1,
879  $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
880  $ iinfo )
881  IF( iinfo.NE.0 ) THEN
882  result( 16 ) = ulpinv
883  result( 17 ) = ulpinv
884  WRITE( nounit, fmt = 9999 )'CGEESX9', iinfo, n, iseed( 1 )
885  info = abs( iinfo )
886  GO TO 270
887  END IF
888 *
889 * Compare condition number for average of selected eigenvalues
890 * taking its condition number into account
891 *
892  anorm = clange( '1', n, n, a, lda, rwork )
893  v = max( REAL( n )*EPS*ANORM, SMLNUM )
894  IF( anorm.EQ.zero )
895  $ v = one
896  IF( v.GT.rcondv ) THEN
897  tol = one
898  ELSE
899  tol = v / rcondv
900  END IF
901  IF( v.GT.rcdvin ) THEN
902  tolin = one
903  ELSE
904  tolin = v / rcdvin
905  END IF
906  tol = max( tol, smlnum / eps )
907  tolin = max( tolin, smlnum / eps )
908  IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
909  result( 16 ) = ulpinv
910  ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
911  result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
912  ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
913  result( 16 ) = ulpinv
914  ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
915  result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
916  ELSE
917  result( 16 ) = one
918  END IF
919 *
920 * Compare condition numbers for right invariant subspace
921 * taking its condition number into account
922 *
923  IF( v.GT.rcondv*rconde ) THEN
924  tol = rcondv
925  ELSE
926  tol = v / rconde
927  END IF
928  IF( v.GT.rcdvin*rcdein ) THEN
929  tolin = rcdvin
930  ELSE
931  tolin = v / rcdein
932  END IF
933  tol = max( tol, smlnum / eps )
934  tolin = max( tolin, smlnum / eps )
935  IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
936  result( 17 ) = ulpinv
937  ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
938  result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
939  ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
940  result( 17 ) = ulpinv
941  ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
942  result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
943  ELSE
944  result( 17 ) = one
945  END IF
946 *
947  270 CONTINUE
948 *
949  END IF
950 *
951  9999 FORMAT( ' CGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
952  $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
953  9998 FORMAT( ' CGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
954  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
955 *
956  RETURN
957 *
958 * End of CGET24
959 *
960  END
cget24
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:337
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cgeesx
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:241
clacpy
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
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cunt01
subroutine cunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
CUNT01
Definition: cunt01.f:128