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