LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
subroutine sgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: sgeesx.f:283
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:118
subroutine sget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, IWORK, BWORK, INFO)
SGET24
Definition: sget24.f:345
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53