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