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