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