LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cget23.f
Go to the documentation of this file.
1*> \brief \b CGET23
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 CGET23( COMP, ISRT, BALANC, JTYPE, THRESH, ISEED,
12* NOUNIT, N, A, LDA, H, W, W1, VL, LDVL, VR,
13* LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
14* RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
15* WORK, LWORK, RWORK, INFO )
16*
17* .. Scalar Arguments ..
18* LOGICAL COMP
19* CHARACTER BALANC
20* INTEGER INFO, ISRT, JTYPE, LDA, LDLRE, LDVL, LDVR,
21* $ LWORK, N, NOUNIT
22* REAL THRESH
23* ..
24* .. Array Arguments ..
25* INTEGER ISEED( 4 )
26* REAL RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
27* $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
28* $ RESULT( 11 ), RWORK( * ), SCALE( * ),
29* $ SCALE1( * )
30* COMPLEX A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
31* $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ),
32* $ WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> CGET23 checks the nonsymmetric eigenvalue problem driver CGEEVX.
42*> If COMP = .FALSE., the first 8 of the following tests will be
43*> performed on the input matrix A, and also test 9 if LWORK is
44*> sufficiently large.
45*> if COMP is .TRUE. all 11 tests will be performed.
46*>
47*> (1) | A * VR - VR * W | / ( n |A| ulp )
48*>
49*> Here VR is the matrix of unit right eigenvectors.
50*> W is a diagonal matrix with diagonal entries W(j).
51*>
52*> (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
53*>
54*> Here VL is the matrix of unit left eigenvectors, A**H is the
55*> conjugate transpose of A, and W is as above.
56*>
57*> (3) | |VR(i)| - 1 | / ulp and largest component real
58*>
59*> VR(i) denotes the i-th column of VR.
60*>
61*> (4) | |VL(i)| - 1 | / ulp and largest component real
62*>
63*> VL(i) denotes the i-th column of VL.
64*>
65*> (5) 0 if W(full) = W(partial), 1/ulp otherwise
66*>
67*> W(full) denotes the eigenvalues computed when VR, VL, RCONDV
68*> and RCONDE are also computed, and W(partial) denotes the
69*> eigenvalues computed when only some of VR, VL, RCONDV, and
70*> RCONDE are computed.
71*>
72*> (6) 0 if VR(full) = VR(partial), 1/ulp otherwise
73*>
74*> VR(full) denotes the right eigenvectors computed when VL, RCONDV
75*> and RCONDE are computed, and VR(partial) denotes the result
76*> when only some of VL and RCONDV are computed.
77*>
78*> (7) 0 if VL(full) = VL(partial), 1/ulp otherwise
79*>
80*> VL(full) denotes the left eigenvectors computed when VR, RCONDV
81*> and RCONDE are computed, and VL(partial) denotes the result
82*> when only some of VR and RCONDV are computed.
83*>
84*> (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
85*> SCALE, ILO, IHI, ABNRM (partial)
86*> 1/ulp otherwise
87*>
88*> SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
89*> (full) is when VR, VL, RCONDE and RCONDV are also computed, and
90*> (partial) is when some are not computed.
91*>
92*> (9) 0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise
93*>
94*> RCONDV(full) denotes the reciprocal condition numbers of the
95*> right eigenvectors computed when VR, VL and RCONDE are also
96*> computed. RCONDV(partial) denotes the reciprocal condition
97*> numbers when only some of VR, VL and RCONDE are computed.
98*>
99*> (10) |RCONDV - RCDVIN| / cond(RCONDV)
100*>
101*> RCONDV is the reciprocal right eigenvector condition number
102*> computed by CGEEVX and RCDVIN (the precomputed true value)
103*> is supplied as input. cond(RCONDV) is the condition number of
104*> RCONDV, and takes errors in computing RCONDV into account, so
105*> that the resulting quantity should be O(ULP). cond(RCONDV) is
106*> essentially given by norm(A)/RCONDE.
107*>
108*> (11) |RCONDE - RCDEIN| / cond(RCONDE)
109*>
110*> RCONDE is the reciprocal eigenvalue condition number
111*> computed by CGEEVX and RCDEIN (the precomputed true value)
112*> is supplied as input. cond(RCONDE) is the condition number
113*> of RCONDE, and takes errors in computing RCONDE into account,
114*> so that the resulting quantity should be O(ULP). cond(RCONDE)
115*> is essentially given by norm(A)/RCONDV.
116*> \endverbatim
117*
118* Arguments:
119* ==========
120*
121*> \param[in] COMP
122*> \verbatim
123*> COMP is LOGICAL
124*> COMP describes which input tests to perform:
125*> = .FALSE. if the computed condition numbers are not to
126*> be tested against RCDVIN and RCDEIN
127*> = .TRUE. if they are to be compared
128*> \endverbatim
129*>
130*> \param[in] ISRT
131*> \verbatim
132*> ISRT is INTEGER
133*> If COMP = .TRUE., ISRT indicates in how the eigenvalues
134*> corresponding to values in RCDVIN and RCDEIN are ordered:
135*> = 0 means the eigenvalues are sorted by
136*> increasing real part
137*> = 1 means the eigenvalues are sorted by
138*> increasing imaginary part
139*> If COMP = .FALSE., ISRT is not referenced.
140*> \endverbatim
141*>
142*> \param[in] BALANC
143*> \verbatim
144*> BALANC is CHARACTER
145*> Describes the balancing option to be tested.
146*> = 'N' for no permuting or diagonal scaling
147*> = 'P' for permuting but no diagonal scaling
148*> = 'S' for no permuting but diagonal scaling
149*> = 'B' for permuting and diagonal scaling
150*> \endverbatim
151*>
152*> \param[in] JTYPE
153*> \verbatim
154*> JTYPE is INTEGER
155*> Type of input matrix. Used to label output if error occurs.
156*> \endverbatim
157*>
158*> \param[in] THRESH
159*> \verbatim
160*> THRESH is REAL
161*> A test will count as "failed" if the "error", computed as
162*> described above, exceeds THRESH. Note that the error
163*> is scaled to be O(1), so THRESH should be a reasonably
164*> small multiple of 1, e.g., 10 or 100. In particular,
165*> it should not depend on the precision (single vs. double)
166*> or the size of the matrix. It must be at least zero.
167*> \endverbatim
168*>
169*> \param[in] ISEED
170*> \verbatim
171*> ISEED is INTEGER array, dimension (4)
172*> If COMP = .FALSE., the random number generator seed
173*> used to produce matrix.
174*> If COMP = .TRUE., ISEED(1) = the number of the example.
175*> Used to label output if error occurs.
176*> \endverbatim
177*>
178*> \param[in] NOUNIT
179*> \verbatim
180*> NOUNIT is INTEGER
181*> The FORTRAN unit number for printing out error messages
182*> (e.g., if a routine returns INFO not equal to 0.)
183*> \endverbatim
184*>
185*> \param[in] N
186*> \verbatim
187*> N is INTEGER
188*> The dimension of A. N must be at least 0.
189*> \endverbatim
190*>
191*> \param[in,out] A
192*> \verbatim
193*> A is COMPLEX array, dimension (LDA,N)
194*> Used to hold the matrix whose eigenvalues are to be
195*> computed.
196*> \endverbatim
197*>
198*> \param[in] LDA
199*> \verbatim
200*> LDA is INTEGER
201*> The leading dimension of A, and H. LDA must be at
202*> least 1 and at least N.
203*> \endverbatim
204*>
205*> \param[out] H
206*> \verbatim
207*> H is COMPLEX array, dimension (LDA,N)
208*> Another copy of the test matrix A, modified by CGEEVX.
209*> \endverbatim
210*>
211*> \param[out] W
212*> \verbatim
213*> W is COMPLEX array, dimension (N)
214*> Contains the eigenvalues of A.
215*> \endverbatim
216*>
217*> \param[out] W1
218*> \verbatim
219*> W1 is COMPLEX array, dimension (N)
220*> Like W, this array contains the eigenvalues of A,
221*> but those computed when CGEEVX only computes a partial
222*> eigendecomposition, i.e. not the eigenvalues and left
223*> and right eigenvectors.
224*> \endverbatim
225*>
226*> \param[out] VL
227*> \verbatim
228*> VL is COMPLEX array, dimension (LDVL,N)
229*> VL holds the computed left eigenvectors.
230*> \endverbatim
231*>
232*> \param[in] LDVL
233*> \verbatim
234*> LDVL is INTEGER
235*> Leading dimension of VL. Must be at least max(1,N).
236*> \endverbatim
237*>
238*> \param[out] VR
239*> \verbatim
240*> VR is COMPLEX array, dimension (LDVR,N)
241*> VR holds the computed right eigenvectors.
242*> \endverbatim
243*>
244*> \param[in] LDVR
245*> \verbatim
246*> LDVR is INTEGER
247*> Leading dimension of VR. Must be at least max(1,N).
248*> \endverbatim
249*>
250*> \param[out] LRE
251*> \verbatim
252*> LRE is COMPLEX array, dimension (LDLRE,N)
253*> LRE holds the computed right or left eigenvectors.
254*> \endverbatim
255*>
256*> \param[in] LDLRE
257*> \verbatim
258*> LDLRE is INTEGER
259*> Leading dimension of LRE. Must be at least max(1,N).
260*> \endverbatim
261*>
262*> \param[out] RCONDV
263*> \verbatim
264*> RCONDV is REAL array, dimension (N)
265*> RCONDV holds the computed reciprocal condition numbers
266*> for eigenvectors.
267*> \endverbatim
268*>
269*> \param[out] RCNDV1
270*> \verbatim
271*> RCNDV1 is REAL array, dimension (N)
272*> RCNDV1 holds more computed reciprocal condition numbers
273*> for eigenvectors.
274*> \endverbatim
275*>
276*> \param[in] RCDVIN
277*> \verbatim
278*> RCDVIN is REAL array, dimension (N)
279*> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
280*> condition numbers for eigenvectors to be compared with
281*> RCONDV.
282*> \endverbatim
283*>
284*> \param[out] RCONDE
285*> \verbatim
286*> RCONDE is REAL array, dimension (N)
287*> RCONDE holds the computed reciprocal condition numbers
288*> for eigenvalues.
289*> \endverbatim
290*>
291*> \param[out] RCNDE1
292*> \verbatim
293*> RCNDE1 is REAL array, dimension (N)
294*> RCNDE1 holds more computed reciprocal condition numbers
295*> for eigenvalues.
296*> \endverbatim
297*>
298*> \param[in] RCDEIN
299*> \verbatim
300*> RCDEIN is REAL array, dimension (N)
301*> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
302*> condition numbers for eigenvalues to be compared with
303*> RCONDE.
304*> \endverbatim
305*>
306*> \param[out] SCALE
307*> \verbatim
308*> SCALE is REAL array, dimension (N)
309*> Holds information describing balancing of matrix.
310*> \endverbatim
311*>
312*> \param[out] SCALE1
313*> \verbatim
314*> SCALE1 is REAL array, dimension (N)
315*> Holds information describing balancing of matrix.
316*> \endverbatim
317*>
318*> \param[out] RESULT
319*> \verbatim
320*> RESULT is REAL array, dimension (11)
321*> The values computed by the 11 tests described above.
322*> The values are currently limited to 1/ulp, to avoid
323*> overflow.
324*> \endverbatim
325*>
326*> \param[out] WORK
327*> \verbatim
328*> WORK is COMPLEX array, dimension (LWORK)
329*> \endverbatim
330*>
331*> \param[in] LWORK
332*> \verbatim
333*> LWORK is INTEGER
334*> The number of entries in WORK. This must be at least
335*> 2*N, and 2*N+N**2 if tests 9, 10 or 11 are to be performed.
336*> \endverbatim
337*>
338*> \param[out] RWORK
339*> \verbatim
340*> RWORK is REAL array, dimension (2*N)
341*> \endverbatim
342*>
343*> \param[out] INFO
344*> \verbatim
345*> INFO is INTEGER
346*> If 0, successful exit.
347*> If <0, input parameter -INFO had an incorrect value.
348*> If >0, CGEEVX returned an error code, the absolute
349*> value of which is returned.
350*> \endverbatim
351*
352* Authors:
353* ========
354*
355*> \author Univ. of Tennessee
356*> \author Univ. of California Berkeley
357*> \author Univ. of Colorado Denver
358*> \author NAG Ltd.
359*
360*> \ingroup complex_eig
361*
362* =====================================================================
363 SUBROUTINE cget23( COMP, ISRT, BALANC, JTYPE, THRESH, ISEED,
364 $ NOUNIT, N, A, LDA, H, W, W1, VL, LDVL, VR,
365 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
366 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
367 $ WORK, LWORK, RWORK, INFO )
368*
369* -- LAPACK test routine --
370* -- LAPACK is a software package provided by Univ. of Tennessee, --
371* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
372*
373* .. Scalar Arguments ..
374 LOGICAL COMP
375 CHARACTER BALANC
376 INTEGER INFO, ISRT, JTYPE, LDA, LDLRE, LDVL, LDVR,
377 $ LWORK, N, NOUNIT
378 REAL THRESH
379* ..
380* .. Array Arguments ..
381 INTEGER ISEED( 4 )
382 REAL RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
383 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
384 $ RESULT( 11 ), RWORK( * ), SCALE( * ),
385 $ scale1( * )
386 COMPLEX A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
387 $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ),
388 $ WORK( * )
389* ..
390*
391* =====================================================================
392*
393* .. Parameters ..
394 REAL ZERO, ONE, TWO
395 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0 )
396 REAL EPSIN
397 PARAMETER ( EPSIN = 5.9605e-8 )
398* ..
399* .. Local Scalars ..
400 LOGICAL BALOK, NOBAL
401 CHARACTER SENSE
402 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
403 $ J, JJ, KMIN
404 REAL ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
405 $ ulp, ulpinv, v, vmax, vmx, vricmp, vrimin,
406 $ vrmx, vtst
407 COMPLEX CTMP
408* ..
409* .. Local Arrays ..
410 CHARACTER SENS( 2 )
411 REAL RES( 2 )
412 COMPLEX CDUM( 1 )
413* ..
414* .. External Functions ..
415 LOGICAL LSAME
416 REAL SCNRM2, SLAMCH
417 EXTERNAL LSAME, SCNRM2, SLAMCH
418* ..
419* .. External Subroutines ..
420 EXTERNAL cgeevx, cget22, clacpy, xerbla
421* ..
422* .. Intrinsic Functions ..
423 INTRINSIC abs, aimag, max, min, real
424* ..
425* .. Data statements ..
426 DATA sens / 'N', 'V' /
427* ..
428* .. Executable Statements ..
429*
430* Check for errors
431*
432 nobal = lsame( balanc, 'N' )
433 balok = nobal .OR. lsame( balanc, 'P' ) .OR.
434 $ lsame( balanc, 'S' ) .OR. lsame( balanc, 'B' )
435 info = 0
436 IF( isrt.NE.0 .AND. isrt.NE.1 ) THEN
437 info = -2
438 ELSE IF( .NOT.balok ) THEN
439 info = -3
440 ELSE IF( thresh.LT.zero ) THEN
441 info = -5
442 ELSE IF( nounit.LE.0 ) THEN
443 info = -7
444 ELSE IF( n.LT.0 ) THEN
445 info = -8
446 ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
447 info = -10
448 ELSE IF( ldvl.LT.1 .OR. ldvl.LT.n ) THEN
449 info = -15
450 ELSE IF( ldvr.LT.1 .OR. ldvr.LT.n ) THEN
451 info = -17
452 ELSE IF( ldlre.LT.1 .OR. ldlre.LT.n ) THEN
453 info = -19
454 ELSE IF( lwork.LT.2*n .OR. ( comp .AND. lwork.LT.2*n+n*n ) ) THEN
455 info = -30
456 END IF
457*
458 IF( info.NE.0 ) THEN
459 CALL xerbla( 'CGET23', -info )
460 RETURN
461 END IF
462*
463* Quick return if nothing to do
464*
465 DO 10 i = 1, 11
466 result( i ) = -one
467 10 CONTINUE
468*
469 IF( n.EQ.0 )
470 $ RETURN
471*
472* More Important constants
473*
474 ulp = slamch( 'Precision' )
475 smlnum = slamch( 'S' )
476 ulpinv = one / ulp
477*
478* Compute eigenvalues and eigenvectors, and test them
479*
480 IF( lwork.GE.2*n+n*n ) THEN
481 sense = 'B'
482 isensm = 2
483 ELSE
484 sense = 'E'
485 isensm = 1
486 END IF
487 CALL clacpy( 'F', n, n, a, lda, h, lda )
488 CALL cgeevx( balanc, 'V', 'V', sense, n, h, lda, w, vl, ldvl, vr,
489 $ ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work,
490 $ lwork, rwork, iinfo )
491 IF( iinfo.NE.0 ) THEN
492 result( 1 ) = ulpinv
493 IF( jtype.NE.22 ) THEN
494 WRITE( nounit, fmt = 9998 )'CGEEVX1', iinfo, n, jtype,
495 $ balanc, iseed
496 ELSE
497 WRITE( nounit, fmt = 9999 )'CGEEVX1', iinfo, n, iseed( 1 )
498 END IF
499 info = abs( iinfo )
500 RETURN
501 END IF
502*
503* Do Test (1)
504*
505 CALL cget22( 'N', 'N', 'N', n, a, lda, vr, ldvr, w, work, rwork,
506 $ res )
507 result( 1 ) = res( 1 )
508*
509* Do Test (2)
510*
511 CALL cget22( 'C', 'N', 'C', n, a, lda, vl, ldvl, w, work, rwork,
512 $ res )
513 result( 2 ) = res( 1 )
514*
515* Do Test (3)
516*
517 DO 30 j = 1, n
518 tnrm = scnrm2( n, vr( 1, j ), 1 )
519 result( 3 ) = max( result( 3 ),
520 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
521 vmx = zero
522 vrmx = zero
523 DO 20 jj = 1, n
524 vtst = abs( vr( jj, j ) )
525 IF( vtst.GT.vmx )
526 $ vmx = vtst
527 IF( aimag( vr( jj, j ) ).EQ.zero .AND.
528 $ abs( real( vr( jj, j ) ) ).GT.vrmx )
529 $ vrmx = abs( real( vr( jj, j ) ) )
530 20 CONTINUE
531 IF( vrmx / vmx.LT.one-two*ulp )
532 $ result( 3 ) = ulpinv
533 30 CONTINUE
534*
535* Do Test (4)
536*
537 DO 50 j = 1, n
538 tnrm = scnrm2( n, vl( 1, j ), 1 )
539 result( 4 ) = max( result( 4 ),
540 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
541 vmx = zero
542 vrmx = zero
543 DO 40 jj = 1, n
544 vtst = abs( vl( jj, j ) )
545 IF( vtst.GT.vmx )
546 $ vmx = vtst
547 IF( aimag( vl( jj, j ) ).EQ.zero .AND.
548 $ abs( real( vl( jj, j ) ) ).GT.vrmx )
549 $ vrmx = abs( real( vl( jj, j ) ) )
550 40 CONTINUE
551 IF( vrmx / vmx.LT.one-two*ulp )
552 $ result( 4 ) = ulpinv
553 50 CONTINUE
554*
555* Test for all options of computing condition numbers
556*
557 DO 200 isens = 1, isensm
558*
559 sense = sens( isens )
560*
561* Compute eigenvalues only, and test them
562*
563 CALL clacpy( 'F', n, n, a, lda, h, lda )
564 CALL cgeevx( balanc, 'N', 'N', sense, n, h, lda, w1, cdum, 1,
565 $ cdum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
566 $ rcndv1, work, lwork, rwork, iinfo )
567 IF( iinfo.NE.0 ) THEN
568 result( 1 ) = ulpinv
569 IF( jtype.NE.22 ) THEN
570 WRITE( nounit, fmt = 9998 )'CGEEVX2', iinfo, n, jtype,
571 $ balanc, iseed
572 ELSE
573 WRITE( nounit, fmt = 9999 )'CGEEVX2', iinfo, n,
574 $ iseed( 1 )
575 END IF
576 info = abs( iinfo )
577 GO TO 190
578 END IF
579*
580* Do Test (5)
581*
582 DO 60 j = 1, n
583 IF( w( j ).NE.w1( j ) )
584 $ result( 5 ) = ulpinv
585 60 CONTINUE
586*
587* Do Test (8)
588*
589 IF( .NOT.nobal ) THEN
590 DO 70 j = 1, n
591 IF( scale( j ).NE.scale1( j ) )
592 $ result( 8 ) = ulpinv
593 70 CONTINUE
594 IF( ilo.NE.ilo1 )
595 $ result( 8 ) = ulpinv
596 IF( ihi.NE.ihi1 )
597 $ result( 8 ) = ulpinv
598 IF( abnrm.NE.abnrm1 )
599 $ result( 8 ) = ulpinv
600 END IF
601*
602* Do Test (9)
603*
604 IF( isens.EQ.2 .AND. n.GT.1 ) THEN
605 DO 80 j = 1, n
606 IF( rcondv( j ).NE.rcndv1( j ) )
607 $ result( 9 ) = ulpinv
608 80 CONTINUE
609 END IF
610*
611* Compute eigenvalues and right eigenvectors, and test them
612*
613 CALL clacpy( 'F', n, n, a, lda, h, lda )
614 CALL cgeevx( balanc, 'N', 'V', sense, n, h, lda, w1, cdum, 1,
615 $ lre, ldlre, ilo1, ihi1, scale1, abnrm1, rcnde1,
616 $ rcndv1, work, lwork, rwork, iinfo )
617 IF( iinfo.NE.0 ) THEN
618 result( 1 ) = ulpinv
619 IF( jtype.NE.22 ) THEN
620 WRITE( nounit, fmt = 9998 )'CGEEVX3', iinfo, n, jtype,
621 $ balanc, iseed
622 ELSE
623 WRITE( nounit, fmt = 9999 )'CGEEVX3', iinfo, n,
624 $ iseed( 1 )
625 END IF
626 info = abs( iinfo )
627 GO TO 190
628 END IF
629*
630* Do Test (5) again
631*
632 DO 90 j = 1, n
633 IF( w( j ).NE.w1( j ) )
634 $ result( 5 ) = ulpinv
635 90 CONTINUE
636*
637* Do Test (6)
638*
639 DO 110 j = 1, n
640 DO 100 jj = 1, n
641 IF( vr( j, jj ).NE.lre( j, jj ) )
642 $ result( 6 ) = ulpinv
643 100 CONTINUE
644 110 CONTINUE
645*
646* Do Test (8) again
647*
648 IF( .NOT.nobal ) THEN
649 DO 120 j = 1, n
650 IF( scale( j ).NE.scale1( j ) )
651 $ result( 8 ) = ulpinv
652 120 CONTINUE
653 IF( ilo.NE.ilo1 )
654 $ result( 8 ) = ulpinv
655 IF( ihi.NE.ihi1 )
656 $ result( 8 ) = ulpinv
657 IF( abnrm.NE.abnrm1 )
658 $ result( 8 ) = ulpinv
659 END IF
660*
661* Do Test (9) again
662*
663 IF( isens.EQ.2 .AND. n.GT.1 ) THEN
664 DO 130 j = 1, n
665 IF( rcondv( j ).NE.rcndv1( j ) )
666 $ result( 9 ) = ulpinv
667 130 CONTINUE
668 END IF
669*
670* Compute eigenvalues and left eigenvectors, and test them
671*
672 CALL clacpy( 'F', n, n, a, lda, h, lda )
673 CALL cgeevx( balanc, 'V', 'N', sense, n, h, lda, w1, lre,
674 $ ldlre, cdum, 1, ilo1, ihi1, scale1, abnrm1,
675 $ rcnde1, rcndv1, work, lwork, rwork, iinfo )
676 IF( iinfo.NE.0 ) THEN
677 result( 1 ) = ulpinv
678 IF( jtype.NE.22 ) THEN
679 WRITE( nounit, fmt = 9998 )'CGEEVX4', iinfo, n, jtype,
680 $ balanc, iseed
681 ELSE
682 WRITE( nounit, fmt = 9999 )'CGEEVX4', iinfo, n,
683 $ iseed( 1 )
684 END IF
685 info = abs( iinfo )
686 GO TO 190
687 END IF
688*
689* Do Test (5) again
690*
691 DO 140 j = 1, n
692 IF( w( j ).NE.w1( j ) )
693 $ result( 5 ) = ulpinv
694 140 CONTINUE
695*
696* Do Test (7)
697*
698 DO 160 j = 1, n
699 DO 150 jj = 1, n
700 IF( vl( j, jj ).NE.lre( j, jj ) )
701 $ result( 7 ) = ulpinv
702 150 CONTINUE
703 160 CONTINUE
704*
705* Do Test (8) again
706*
707 IF( .NOT.nobal ) THEN
708 DO 170 j = 1, n
709 IF( scale( j ).NE.scale1( j ) )
710 $ result( 8 ) = ulpinv
711 170 CONTINUE
712 IF( ilo.NE.ilo1 )
713 $ result( 8 ) = ulpinv
714 IF( ihi.NE.ihi1 )
715 $ result( 8 ) = ulpinv
716 IF( abnrm.NE.abnrm1 )
717 $ result( 8 ) = ulpinv
718 END IF
719*
720* Do Test (9) again
721*
722 IF( isens.EQ.2 .AND. n.GT.1 ) THEN
723 DO 180 j = 1, n
724 IF( rcondv( j ).NE.rcndv1( j ) )
725 $ result( 9 ) = ulpinv
726 180 CONTINUE
727 END IF
728*
729 190 CONTINUE
730*
731 200 CONTINUE
732*
733* If COMP, compare condition numbers to precomputed ones
734*
735 IF( comp ) THEN
736 CALL clacpy( 'F', n, n, a, lda, h, lda )
737 CALL cgeevx( 'N', 'V', 'V', 'B', n, h, lda, w, vl, ldvl, vr,
738 $ ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
739 $ work, lwork, rwork, iinfo )
740 IF( iinfo.NE.0 ) THEN
741 result( 1 ) = ulpinv
742 WRITE( nounit, fmt = 9999 )'CGEEVX5', iinfo, n, iseed( 1 )
743 info = abs( iinfo )
744 GO TO 250
745 END IF
746*
747* Sort eigenvalues and condition numbers lexicographically
748* to compare with inputs
749*
750 DO 220 i = 1, n - 1
751 kmin = i
752 IF( isrt.EQ.0 ) THEN
753 vrimin = real( w( i ) )
754 ELSE
755 vrimin = aimag( w( i ) )
756 END IF
757 DO 210 j = i + 1, n
758 IF( isrt.EQ.0 ) THEN
759 vricmp = real( w( j ) )
760 ELSE
761 vricmp = aimag( w( j ) )
762 END IF
763 IF( vricmp.LT.vrimin ) THEN
764 kmin = j
765 vrimin = vricmp
766 END IF
767 210 CONTINUE
768 ctmp = w( kmin )
769 w( kmin ) = w( i )
770 w( i ) = ctmp
771 vrimin = rconde( kmin )
772 rconde( kmin ) = rconde( i )
773 rconde( i ) = vrimin
774 vrimin = rcondv( kmin )
775 rcondv( kmin ) = rcondv( i )
776 rcondv( i ) = vrimin
777 220 CONTINUE
778*
779* Compare condition numbers for eigenvectors
780* taking their condition numbers into account
781*
782 result( 10 ) = zero
783 eps = max( epsin, ulp )
784 v = max( real( n )*eps*abnrm, smlnum )
785 IF( abnrm.EQ.zero )
786 $ v = one
787 DO 230 i = 1, n
788 IF( v.GT.rcondv( i )*rconde( i ) ) THEN
789 tol = rcondv( i )
790 ELSE
791 tol = v / rconde( i )
792 END IF
793 IF( v.GT.rcdvin( i )*rcdein( i ) ) THEN
794 tolin = rcdvin( i )
795 ELSE
796 tolin = v / rcdein( i )
797 END IF
798 tol = max( tol, smlnum / eps )
799 tolin = max( tolin, smlnum / eps )
800 IF( eps*( rcdvin( i )-tolin ).GT.rcondv( i )+tol ) THEN
801 vmax = one / eps
802 ELSE IF( rcdvin( i )-tolin.GT.rcondv( i )+tol ) THEN
803 vmax = ( rcdvin( i )-tolin ) / ( rcondv( i )+tol )
804 ELSE IF( rcdvin( i )+tolin.LT.eps*( rcondv( i )-tol ) ) THEN
805 vmax = one / eps
806 ELSE IF( rcdvin( i )+tolin.LT.rcondv( i )-tol ) THEN
807 vmax = ( rcondv( i )-tol ) / ( rcdvin( i )+tolin )
808 ELSE
809 vmax = one
810 END IF
811 result( 10 ) = max( result( 10 ), vmax )
812 230 CONTINUE
813*
814* Compare condition numbers for eigenvalues
815* taking their condition numbers into account
816*
817 result( 11 ) = zero
818 DO 240 i = 1, n
819 IF( v.GT.rcondv( i ) ) THEN
820 tol = one
821 ELSE
822 tol = v / rcondv( i )
823 END IF
824 IF( v.GT.rcdvin( i ) ) THEN
825 tolin = one
826 ELSE
827 tolin = v / rcdvin( i )
828 END IF
829 tol = max( tol, smlnum / eps )
830 tolin = max( tolin, smlnum / eps )
831 IF( eps*( rcdein( i )-tolin ).GT.rconde( i )+tol ) THEN
832 vmax = one / eps
833 ELSE IF( rcdein( i )-tolin.GT.rconde( i )+tol ) THEN
834 vmax = ( rcdein( i )-tolin ) / ( rconde( i )+tol )
835 ELSE IF( rcdein( i )+tolin.LT.eps*( rconde( i )-tol ) ) THEN
836 vmax = one / eps
837 ELSE IF( rcdein( i )+tolin.LT.rconde( i )-tol ) THEN
838 vmax = ( rconde( i )-tol ) / ( rcdein( i )+tolin )
839 ELSE
840 vmax = one
841 END IF
842 result( 11 ) = max( result( 11 ), vmax )
843 240 CONTINUE
844 250 CONTINUE
845*
846 END IF
847*
848 9999 FORMAT( ' CGET23: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
849 $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
850 9998 FORMAT( ' CGET23: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
851 $ i6, ', JTYPE=', i6, ', BALANC = ', a, ', ISEED=(',
852 $ 3( i5, ',' ), i5, ')' )
853*
854 RETURN
855*
856* End of CGET23
857*
858 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cget23(COMP, ISRT, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, W, W1, VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT, WORK, LWORK, RWORK, INFO)
CGET23
Definition: cget23.f:368
subroutine cget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, WORK, RWORK, RESULT)
CGET22
Definition: cget22.f:144
subroutine cgeevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV, WORK, LWORK, RWORK, INFO)
CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: cgeevx.f:288
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