365 SUBROUTINE cget23( COMP, ISRT, BALANC, JTYPE, THRESH, ISEED,
366 $ NOUNIT, N, A, LDA, H, W, W1, VL, LDVL, VR,
367 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
368 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
369 $ WORK, LWORK, RWORK, INFO )
379 INTEGER INFO, ISRT, JTYPE, LDA, LDLRE, LDVL, LDVR,
385 REAL RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
386 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
387 $ RESULT( 11 ), RWORK( * ), SCALE( * ),
389 COMPLEX A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
390 $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ),
398 PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0 )
400 PARAMETER ( EPSIN = 5.9605e-8 )
405 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
407 REAL ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
408 $ ulp, ulpinv, v, vmax, vmx, vricmp, vrimin,
420 EXTERNAL LSAME, SCNRM2, SLAMCH
426 INTRINSIC abs, aimag, max, min, real
429 DATA sens /
'N',
'V' /
435 nobal = lsame( balanc,
'N' )
436 balok = nobal .OR. lsame( balanc,
'P' ) .OR.
437 $ lsame( balanc,
'S' ) .OR. lsame( balanc,
'B' )
439 IF( isrt.NE.0 .AND. isrt.NE.1 )
THEN
441 ELSE IF( .NOT.balok )
THEN
443 ELSE IF( thresh.LT.zero )
THEN
445 ELSE IF( nounit.LE.0 )
THEN
447 ELSE IF( n.LT.0 )
THEN
449 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
451 ELSE IF( ldvl.LT.1 .OR. ldvl.LT.n )
THEN
453 ELSE IF( ldvr.LT.1 .OR. ldvr.LT.n )
THEN
455 ELSE IF( ldlre.LT.1 .OR. ldlre.LT.n )
THEN
457 ELSE IF( lwork.LT.2*n .OR. ( comp .AND. lwork.LT.2*n+n*n ) )
THEN
462 CALL xerbla(
'CGET23', -info )
477 ulp = slamch(
'Precision' )
478 smlnum = slamch(
'S' )
483 IF( lwork.GE.2*n+n*n )
THEN
490 CALL clacpy(
'F', n, n, a, lda, h, lda )
491 CALL cgeevx( balanc,
'V',
'V', sense, n, h, lda, w, vl, ldvl, vr,
492 $ ldvr, ilo, ihi, scale, abnrm, rconde, rcondv, work,
493 $ lwork, rwork, iinfo )
494 IF( iinfo.NE.0 )
THEN
496 IF( jtype.NE.22 )
THEN
497 WRITE( nounit, fmt = 9998 )
'CGEEVX1', iinfo, n, jtype,
500 WRITE( nounit, fmt = 9999 )
'CGEEVX1', iinfo, n, iseed( 1 )
508 CALL cget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, w, work, rwork,
510 result( 1 ) = res( 1 )
514 CALL cget22(
'C',
'N',
'C', n, a, lda, vl, ldvl, w, work, rwork,
516 result( 2 ) = res( 1 )
521 tnrm = scnrm2( n, vr( 1, j ), 1 )
522 result( 3 ) = max( result( 3 ),
523 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
527 vtst = abs( vr( jj, j ) )
530 IF( aimag( vr( jj, j ) ).EQ.zero .AND.
531 $ abs( real( vr( jj, j ) ) ).GT.vrmx )
532 $ vrmx = abs( real( vr( jj, j ) ) )
534 IF( vrmx / vmx.LT.one-two*ulp )
535 $ result( 3 ) = ulpinv
541 tnrm = scnrm2( n, vl( 1, j ), 1 )
542 result( 4 ) = max( result( 4 ),
543 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
547 vtst = abs( vl( jj, j ) )
550 IF( aimag( vl( jj, j ) ).EQ.zero .AND.
551 $ abs( real( vl( jj, j ) ) ).GT.vrmx )
552 $ vrmx = abs( real( vl( jj, j ) ) )
554 IF( vrmx / vmx.LT.one-two*ulp )
555 $ result( 4 ) = ulpinv
560 DO 200 isens = 1, isensm
562 sense = sens( isens )
566 CALL clacpy(
'F', n, n, a, lda, h, lda )
567 CALL cgeevx( balanc,
'N',
'N', sense, n, h, lda, w1, cdum, 1,
568 $ cdum, 1, ilo1, ihi1, scale1, abnrm1, rcnde1,
569 $ rcndv1, work, lwork, rwork, iinfo )
570 IF( iinfo.NE.0 )
THEN
572 IF( jtype.NE.22 )
THEN
573 WRITE( nounit, fmt = 9998 )
'CGEEVX2', iinfo, n, jtype,
576 WRITE( nounit, fmt = 9999 )
'CGEEVX2', iinfo, n,
586 IF( w( j ).NE.w1( j ) )
587 $ result( 5 ) = ulpinv
592 IF( .NOT.nobal )
THEN
594 IF( scale( j ).NE.scale1( j ) )
595 $ result( 8 ) = ulpinv
598 $ result( 8 ) = ulpinv
600 $ result( 8 ) = ulpinv
601 IF( abnrm.NE.abnrm1 )
602 $ result( 8 ) = ulpinv
607 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
609 IF( rcondv( j ).NE.rcndv1( j ) )
610 $ result( 9 ) = ulpinv
616 CALL clacpy(
'F', n, n, a, lda, h, lda )
617 CALL cgeevx( balanc,
'N',
'V', sense, n, h, lda, w1, cdum, 1,
618 $ lre, ldlre, ilo1, ihi1, scale1, abnrm1, rcnde1,
619 $ rcndv1, work, lwork, rwork, iinfo )
620 IF( iinfo.NE.0 )
THEN
622 IF( jtype.NE.22 )
THEN
623 WRITE( nounit, fmt = 9998 )
'CGEEVX3', iinfo, n, jtype,
626 WRITE( nounit, fmt = 9999 )
'CGEEVX3', iinfo, n,
636 IF( w( j ).NE.w1( j ) )
637 $ result( 5 ) = ulpinv
644 IF( vr( j, jj ).NE.lre( j, jj ) )
645 $ result( 6 ) = ulpinv
651 IF( .NOT.nobal )
THEN
653 IF( scale( j ).NE.scale1( j ) )
654 $ result( 8 ) = ulpinv
657 $ result( 8 ) = ulpinv
659 $ result( 8 ) = ulpinv
660 IF( abnrm.NE.abnrm1 )
661 $ result( 8 ) = ulpinv
666 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
668 IF( rcondv( j ).NE.rcndv1( j ) )
669 $ result( 9 ) = ulpinv
675 CALL clacpy(
'F', n, n, a, lda, h, lda )
676 CALL cgeevx( balanc,
'V',
'N', sense, n, h, lda, w1, lre,
677 $ ldlre, cdum, 1, ilo1, ihi1, scale1, abnrm1,
678 $ rcnde1, rcndv1, work, lwork, rwork, iinfo )
679 IF( iinfo.NE.0 )
THEN
681 IF( jtype.NE.22 )
THEN
682 WRITE( nounit, fmt = 9998 )
'CGEEVX4', iinfo, n, jtype,
685 WRITE( nounit, fmt = 9999 )
'CGEEVX4', iinfo, n,
695 IF( w( j ).NE.w1( j ) )
696 $ result( 5 ) = ulpinv
703 IF( vl( j, jj ).NE.lre( j, jj ) )
704 $ result( 7 ) = ulpinv
710 IF( .NOT.nobal )
THEN
712 IF( scale( j ).NE.scale1( j ) )
713 $ result( 8 ) = ulpinv
716 $ result( 8 ) = ulpinv
718 $ result( 8 ) = ulpinv
719 IF( abnrm.NE.abnrm1 )
720 $ result( 8 ) = ulpinv
725 IF( isens.EQ.2 .AND. n.GT.1 )
THEN
727 IF( rcondv( j ).NE.rcndv1( j ) )
728 $ result( 9 ) = ulpinv
739 CALL clacpy(
'F', n, n, a, lda, h, lda )
740 CALL cgeevx(
'N',
'V',
'V',
'B', n, h, lda, w, vl, ldvl, vr,
741 $ ldvr, ilo, ihi, scale, abnrm, rconde, rcondv,
742 $ work, lwork, rwork, iinfo )
743 IF( iinfo.NE.0 )
THEN
745 WRITE( nounit, fmt = 9999 )
'CGEEVX5', iinfo, n, iseed( 1 )
756 vrimin = real( w( i ) )
758 vrimin = aimag( w( i ) )
762 vricmp = real( w( j ) )
764 vricmp = aimag( w( j ) )
766 IF( vricmp.LT.vrimin )
THEN
774 vrimin = rconde( kmin )
775 rconde( kmin ) = rconde( i )
777 vrimin = rcondv( kmin )
778 rcondv( kmin ) = rcondv( i )
786 eps = max( epsin, ulp )
787 v = max( real( n )*eps*abnrm, smlnum )
791 IF( v.GT.rcondv( i )*rconde( i ) )
THEN
794 tol = v / rconde( i )
796 IF( v.GT.rcdvin( i )*rcdein( i ) )
THEN
799 tolin = v / rcdein( i )
801 tol = max( tol, smlnum / eps )
802 tolin = max( tolin, smlnum / eps )
803 IF( eps*( rcdvin( i )-tolin ).GT.rcondv( i )+tol )
THEN
805 ELSE IF( rcdvin( i )-tolin.GT.rcondv( i )+tol )
THEN
806 vmax = ( rcdvin( i )-tolin ) / ( rcondv( i )+tol )
807 ELSE IF( rcdvin( i )+tolin.LT.eps*( rcondv( i )-tol ) )
THEN
809 ELSE IF( rcdvin( i )+tolin.LT.rcondv( i )-tol )
THEN
810 vmax = ( rcondv( i )-tol ) / ( rcdvin( i )+tolin )
814 result( 10 ) = max( result( 10 ), vmax )
822 IF( v.GT.rcondv( i ) )
THEN
825 tol = v / rcondv( i )
827 IF( v.GT.rcdvin( i ) )
THEN
830 tolin = v / rcdvin( i )
832 tol = max( tol, smlnum / eps )
833 tolin = max( tolin, smlnum / eps )
834 IF( eps*( rcdein( i )-tolin ).GT.rconde( i )+tol )
THEN
836 ELSE IF( rcdein( i )-tolin.GT.rconde( i )+tol )
THEN
837 vmax = ( rcdein( i )-tolin ) / ( rconde( i )+tol )
838 ELSE IF( rcdein( i )+tolin.LT.eps*( rconde( i )-tol ) )
THEN
840 ELSE IF( rcdein( i )+tolin.LT.rconde( i )-tol )
THEN
841 vmax = ( rconde( i )-tol ) / ( rcdein( i )+tolin )
845 result( 11 ) = max( result( 11 ), vmax )
851 9999
FORMAT(
' CGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
852 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
853 9998
FORMAT(
' CGET23: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
854 $ i6,
', JTYPE=', i6,
', BALANC = ', a,
', ISEED=(',
855 $ 3( i5,
',' ), i5,
')' )