OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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.NE..AND..NE. IF( IINFO0 IINFON+2 ) THEN
830 RESULT( 15 ) = ULPINV
831.NE. IF( JTYPE22 ) 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.NE. IF( RCNDV1RCONDV )
845 $ RESULT( 15 ) = ULPINV
846*
847* Perform tests (10), (11), (12), and (13)
848*
849 DO 240 I = 1, N
850.NE..OR..NE. IF( WR( I )WRT( I ) WI( I )WIT( I ) )
851 $ RESULT( 10 ) = ULPINV
852 DO 230 J = 1, N
853.NE. IF( H( I, J )HT( I, J ) )
854 $ RESULT( 11 ) = ULPINV
855.NE. IF( VS( I, J )VS1( I, J ) )
856 $ RESULT( 12 ) = ULPINV
857 230 CONTINUE
858 240 CONTINUE
859.NE. IF( SDIMSDIM1 )
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.LT. IF( WRTMP( J )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.NE..AND..NE. IF( IINFO0 IINFON+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.EQ. IF( ANORMZERO )
927 $ V = ONE
928.GT. IF( VRCONDV ) THEN
929 TOL = ONE
930 ELSE
931 TOL = V / RCONDV
932 END IF
933.GT. IF( VRCDVIN ) 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.GT. IF( EPS*( RCDEIN-TOLIN )RCONDE+TOL ) THEN
941 RESULT( 16 ) = ULPINV
942.GT. ELSE IF( RCDEIN-TOLINRCONDE+TOL ) THEN
943 RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
944.LT. ELSE IF( RCDEIN+TOLINEPS*( RCONDE-TOL ) ) THEN
945 RESULT( 16 ) = ULPINV
946.LT. ELSE IF( RCDEIN+TOLINRCONDE-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.GT. IF( VRCONDV*RCONDE ) THEN
956 TOL = RCONDV
957 ELSE
958 TOL = V / RCONDE
959 END IF
960.GT. IF( VRCDVIN*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.GT. IF( EPS*( RCDVIN-TOLIN )RCONDV+TOL ) THEN
968 RESULT( 17 ) = ULPINV
969.GT. ELSE IF( RCDVIN-TOLINRCONDV+TOL ) THEN
970 RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
971.LT. ELSE IF( RCDVIN+TOLINEPS*( RCONDV-TOL ) ) THEN
972 RESULT( 17 ) = ULPINV
973.LT. ELSE IF( RCDVIN+TOLINRCONDV-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
slacpy
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
sgeesx
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
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
sort01
subroutine sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01
Definition sort01.f:116
sget24
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21