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 )
342 INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
344 REAL RCDEIN, RCDVIN, THRESH
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( * )
359 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
360 $ cone = ( 1.0e+0, 0.0e+0 ) )
362 parameter( zero = 0.0e+0, one = 1.0e+0 )
364 parameter( epsin = 5.9605e-8 )
368 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
370 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
371 $ smlnum, tol, tolin, ulp, ulpinv, v, vricmp,
381 EXTERNAL CSLECT, CLANGE, SLAMCH
387 INTRINSIC abs, aimag, max, min, real
391 REAL SELWI( 20 ), SELWR( 20 )
394 INTEGER SELDIM, SELOPT
397 COMMON / sslct / selopt, seldim, selval, selwr, selwi
404 IF( thresh.LT.zero )
THEN
406 ELSE IF( nounit.LE.0 )
THEN
408 ELSE IF( n.LT.0 )
THEN
410 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
412 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n )
THEN
414 ELSE IF( lwork.LT.2*n )
THEN
419 CALL xerbla(
'CGET24', -info )
434 smlnum = slamch(
'Safe minimum' )
435 ulp = slamch(
'Precision' )
442 IF( isort.EQ.0 )
THEN
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,
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,
462 WRITE( nounit, fmt = 9999 )
'CGEESX1', iinfo, n,
468 IF( isort.EQ.0 )
THEN
469 CALL ccopy( n, w, 1, wtmp, 1 )
474 result( 1+rsub ) = zero
477 IF( h( i, j ).NE.czero )
478 $ result( 1+rsub ) = ulpinv
486 CALL clacpy(
' ', n, n, a, lda, vs1, ldvs )
490 CALL cgemm(
'No transpose',
'No transpose', n, n, n, cone, vs,
491 $ ldvs, h, lda, czero, ht, lda )
495 CALL cgemm(
'No transpose',
'Conjugate transpose', n, n, n,
496 $ -cone, ht, lda, vs, ldvs, cone, vs1, ldvs )
498 anorm = max( clange(
'1', n, n, a, lda, rwork ), smlnum )
499 wnorm = clange(
'1', n, n, vs1, ldvs, rwork )
501 IF( anorm.GT.wnorm )
THEN
502 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
504 IF( anorm.LT.one )
THEN
505 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
508 result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
515 CALL cunt01(
'Columns', n, n, vs, ldvs, work, lwork, rwork,
520 result( 4+rsub ) = zero
522 IF( h( i, i ).NE.w( i ) )
523 $ result( 4+rsub ) = ulpinv
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,
532 IF( iinfo.NE.0 )
THEN
533 result( 5+rsub ) = ulpinv
534 IF( jtype.NE.22 )
THEN
535 WRITE( nounit, fmt = 9998 )
'CGEESX2', iinfo, n, jtype,
538 WRITE( nounit, fmt = 9999 )
'CGEESX2', iinfo, n,
545 result( 5+rsub ) = zero
548 IF( h( i, j ).NE.ht( i, j ) )
549 $ result( 5+rsub ) = ulpinv
555 result( 6+rsub ) = zero
557 IF( w( i ).NE.wt( i ) )
558 $ result( 6+rsub ) = ulpinv
563 IF( isort.EQ.1 )
THEN
567 IF( cslect( w( i ) ) )
568 $ knteig = knteig + 1
570 IF( cslect( w( i+1 ) ) .AND.
571 $ ( .NOT.cslect( w( i ) ) ) )result( 13 ) = ulpinv
575 $ result( 13 ) = ulpinv
583 IF( lwork.GE.( n*( n+1 ) ) / 2 )
THEN
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,
594 IF( iinfo.NE.0 )
THEN
595 result( 14 ) = ulpinv
596 result( 15 ) = ulpinv
597 IF( jtype.NE.22 )
THEN
598 WRITE( nounit, fmt = 9998 )
'CGEESX3', iinfo, n, jtype,
601 WRITE( nounit, fmt = 9999 )
'CGEESX3', iinfo, n,
611 IF( w( i ).NE.wt( i ) )
612 $ result( 10 ) = ulpinv
614 IF( h( i, j ).NE.ht( i, j ) )
615 $ result( 11 ) = ulpinv
616 IF( vs( i, j ).NE.vs1( i, j ) )
617 $ result( 12 ) = ulpinv
621 $ result( 13 ) = ulpinv
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,
629 IF( iinfo.NE.0 )
THEN
630 result( 14 ) = ulpinv
631 result( 15 ) = ulpinv
632 IF( jtype.NE.22 )
THEN
633 WRITE( nounit, fmt = 9998 )
'CGEESX4', iinfo, n, jtype,
636 WRITE( nounit, fmt = 9999 )
'CGEESX4', iinfo, n,
645 IF( rcnde1.NE.rconde )
646 $ result( 14 ) = ulpinv
647 IF( rcndv1.NE.rcondv )
648 $ result( 15 ) = ulpinv
653 IF( w( i ).NE.wt( i ) )
654 $ result( 10 ) = ulpinv
656 IF( h( i, j ).NE.ht( i, j ) )
657 $ result( 11 ) = ulpinv
658 IF( vs( i, j ).NE.vs1( i, j ) )
659 $ result( 12 ) = ulpinv
663 $ result( 13 ) = ulpinv
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,
671 IF( iinfo.NE.0 )
THEN
672 result( 14 ) = ulpinv
673 IF( jtype.NE.22 )
THEN
674 WRITE( nounit, fmt = 9998 )
'CGEESX5', iinfo, n, jtype,
677 WRITE( nounit, fmt = 9999 )
'CGEESX5', iinfo, n,
686 IF( rcnde1.NE.rconde )
687 $ result( 14 ) = ulpinv
692 IF( w( i ).NE.wt( i ) )
693 $ result( 10 ) = ulpinv
695 IF( h( i, j ).NE.ht( i, j ) )
696 $ result( 11 ) = ulpinv
697 IF( vs( i, j ).NE.vs1( i, j ) )
698 $ result( 12 ) = ulpinv
702 $ result( 13 ) = ulpinv
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,
710 IF( iinfo.NE.0 )
THEN
711 result( 14 ) = ulpinv
712 IF( jtype.NE.22 )
THEN
713 WRITE( nounit, fmt = 9998 )
'CGEESX6', iinfo, n, jtype,
716 WRITE( nounit, fmt = 9999 )
'CGEESX6', iinfo, n,
725 IF( rcnde1.NE.rconde )
726 $ result( 14 ) = ulpinv
731 IF( w( i ).NE.wt( i ) )
732 $ result( 10 ) = ulpinv
734 IF( h( i, j ).NE.ht( i, j ) )
735 $ result( 11 ) = ulpinv
736 IF( vs( i, j ).NE.vs1( i, j ) )
737 $ result( 12 ) = ulpinv
741 $ result( 13 ) = ulpinv
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,
749 IF( iinfo.NE.0 )
THEN
750 result( 15 ) = ulpinv
751 IF( jtype.NE.22 )
THEN
752 WRITE( nounit, fmt = 9998 )
'CGEESX7', iinfo, n, jtype,
755 WRITE( nounit, fmt = 9999 )
'CGEESX7', iinfo, n,
764 IF( rcndv1.NE.rcondv )
765 $ result( 15 ) = ulpinv
770 IF( w( i ).NE.wt( i ) )
771 $ result( 10 ) = ulpinv
773 IF( h( i, j ).NE.ht( i, j ) )
774 $ result( 11 ) = ulpinv
775 IF( vs( i, j ).NE.vs1( i, j ) )
776 $ result( 12 ) = ulpinv
780 $ result( 13 ) = ulpinv
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,
788 IF( iinfo.NE.0 )
THEN
789 result( 15 ) = ulpinv
790 IF( jtype.NE.22 )
THEN
791 WRITE( nounit, fmt = 9998 )
'CGEESX8', iinfo, n, jtype,
794 WRITE( nounit, fmt = 9999 )
'CGEESX8', iinfo, n,
803 IF( rcndv1.NE.rcondv )
804 $ result( 15 ) = ulpinv
809 IF( w( i ).NE.wt( i ) )
810 $ result( 10 ) = ulpinv
812 IF( h( i, j ).NE.ht( i, j ) )
813 $ result( 11 ) = ulpinv
814 IF( vs( i, j ).NE.vs1( i, j ) )
815 $ result( 12 ) = ulpinv
819 $ result( 13 ) = ulpinv
836 eps = max( ulp, epsin )
839 selval( i ) = .false.
840 selwr( i ) = real( wtmp( i ) )
841 selwi( i ) = aimag( wtmp( i ) )
846 vrimin = real( wtmp( i ) )
848 vrimin = aimag( wtmp( i ) )
852 vricmp = real( wtmp( j ) )
854 vricmp = aimag( wtmp( j ) )
856 IF( vricmp.LT.vrimin )
THEN
862 wtmp( kmin ) = wtmp( i )
865 ipnt( i ) = ipnt( kmin )
869 selval( ipnt( islct( i ) ) ) = .true.
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,
878 IF( iinfo.NE.0 )
THEN
879 result( 16 ) = ulpinv
880 result( 17 ) = ulpinv
881 WRITE( nounit, fmt = 9999 )
'CGEESX9', iinfo, n, iseed( 1 )
889 anorm = clange(
'1', n, n, a, lda, rwork )
890 v = max( real( n )*eps*anorm, smlnum )
893 IF( v.GT.rcondv )
THEN
898 IF( v.GT.rcdvin )
THEN
903 tol = max( tol, smlnum / eps )
904 tolin = max( tolin, smlnum / eps )
905 IF( eps*( rcdein-tolin ).GT.rconde+tol )
THEN
906 result( 16 ) = ulpinv
907 ELSE IF( rcdein-tolin.GT.rconde+tol )
THEN
908 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
909 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) )
THEN
910 result( 16 ) = ulpinv
911 ELSE IF( rcdein+tolin.LT.rconde-tol )
THEN
912 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
920 IF( v.GT.rcondv*rconde )
THEN
925 IF( v.GT.rcdvin*rcdein )
THEN
930 tol = max( tol, smlnum / eps )
931 tolin = max( tolin, smlnum / eps )
932 IF( eps*( rcdvin-tolin ).GT.rcondv+tol )
THEN
933 result( 17 ) = ulpinv
934 ELSE IF( rcdvin-tolin.GT.rcondv+tol )
THEN
935 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
936 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) )
THEN
937 result( 17 ) = ulpinv
938 ELSE IF( rcdvin+tolin.LT.rcondv-tol )
THEN
939 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
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,
')' )