341 SUBROUTINE sget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
342 $ h, ht, wr, wi, wrt, wit, wrtmp, witmp, vs,
343 $ ldvs, vs1, rcdein, rcdvin, nslct, islct,
344 $ result, work, lwork, iwork, bwork, info )
353 INTEGER info, jtype, lda, ldvs, lwork, n, nounit, nslct
354 REAL rcdein, rcdvin, thresh
358 INTEGER iseed( 4 ), islct( * ), iwork( * )
359 REAL a( lda, * ), h( lda, * ), ht( lda, * ),
360 $ result( 17 ), vs( ldvs, * ), vs1( ldvs, * ),
361 $ wi( * ), wit( * ), witmp( * ), work( * ),
362 $ wr( * ), wrt( * ), wrtmp( * )
369 parameter( zero = 0.0e0, one = 1.0e0 )
371 parameter( epsin = 5.9605e-8 )
375 INTEGER i, iinfo, isort, itmp, j, kmin, knteig, liwork,
377 REAL anorm, eps, rcnde1, rcndv1, rconde, rcondv,
378 $ smlnum, tmp, tol, tolin, ulp, ulpinv, v, vimin,
386 REAL selwi( 20 ), selwr( 20 )
389 INTEGER seldim, selopt
392 common / sslct / selopt, seldim, selval, selwr, selwi
403 INTRINSIC abs, max, min,
REAL, sign, sqrt
410 IF( thresh.LT.zero )
THEN
412 ELSE IF( nounit.LE.0 )
THEN
414 ELSE IF( n.LT.0 )
THEN
416 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
418 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n )
THEN
420 ELSE IF( lwork.LT.3*n )
THEN
425 CALL
xerbla(
'SGET24', -info )
440 smlnum =
slamch(
'Safe minimum' )
441 ulp =
slamch(
'Precision' )
449 IF( isort.EQ.0 )
THEN
459 CALL
slacpy(
'F', n, n, a, lda, h, lda )
460 CALL
sgeesx(
'V', sort,
sslect,
'N', n, h, lda, sdim, wr, wi,
461 $ vs, ldvs, rconde, rcondv, work, lwork, iwork,
462 $ liwork, bwork, iinfo )
463 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
464 result( 1+rsub ) = ulpinv
465 IF( jtype.NE.22 )
THEN
466 WRITE( nounit, fmt = 9998 )
'SGEESX1', iinfo, n, jtype,
469 WRITE( nounit, fmt = 9999 )
'SGEESX1', iinfo, n,
475 IF( isort.EQ.0 )
THEN
476 CALL
scopy( n, wr, 1, wrtmp, 1 )
477 CALL
scopy( n, wi, 1, witmp, 1 )
482 result( 1+rsub ) = zero
485 IF( h( i, j ).NE.zero )
486 $ result( 1+rsub ) = ulpinv
490 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
491 $ result( 1+rsub ) = ulpinv
494 IF( h( i+1, i ).NE.zero )
THEN
495 IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
496 $ zero .OR. sign( one, h( i+1, i ) ).EQ.
497 $ sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
505 CALL
slacpy(
' ', n, n, a, lda, vs1, ldvs )
509 CALL
sgemm(
'No transpose',
'No transpose', n, n, n, one, vs,
510 $ ldvs, h, lda, zero, ht, lda )
514 CALL
sgemm(
'No transpose',
'Transpose', n, n, n, -one, ht,
515 $ lda, vs, ldvs, one, vs1, ldvs )
517 anorm = max(
slange(
'1', n, n, a, lda, work ), smlnum )
518 wnorm =
slange(
'1', n, n, vs1, ldvs, work )
520 IF( anorm.GT.wnorm )
THEN
521 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
523 IF( anorm.LT.one )
THEN
524 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
527 result( 2+rsub ) = min( wnorm / anorm,
REAL( N ) ) /
534 CALL
sort01(
'Columns', n, n, vs, ldvs, work, lwork,
539 result( 4+rsub ) = zero
541 IF( h( i, i ).NE.wr( i ) )
542 $ result( 4+rsub ) = ulpinv
545 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
546 $ result( 4+rsub ) = ulpinv
547 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
548 $ result( 4+rsub ) = ulpinv
551 IF( h( i+1, i ).NE.zero )
THEN
552 tmp = sqrt( abs( h( i+1, i ) ) )*
553 $ sqrt( abs( h( i, i+1 ) ) )
554 result( 4+rsub ) = max( result( 4+rsub ),
555 $ abs( wi( i )-tmp ) /
556 $ max( ulp*tmp, smlnum ) )
557 result( 4+rsub ) = max( result( 4+rsub ),
558 $ abs( wi( i+1 )+tmp ) /
559 $ max( ulp*tmp, smlnum ) )
560 ELSE IF( i.GT.1 )
THEN
561 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
562 $ wi( i ).NE.zero )result( 4+rsub ) = ulpinv
568 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
569 CALL
sgeesx(
'N', sort,
sslect,
'N', n, ht, lda, sdim, wrt,
570 $ wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
571 $ liwork, bwork, iinfo )
572 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
573 result( 5+rsub ) = ulpinv
574 IF( jtype.NE.22 )
THEN
575 WRITE( nounit, fmt = 9998 )
'SGEESX2', iinfo, n, jtype,
578 WRITE( nounit, fmt = 9999 )
'SGEESX2', iinfo, n,
585 result( 5+rsub ) = zero
588 IF( h( i, j ).NE.ht( i, j ) )
589 $ result( 5+rsub ) = ulpinv
595 result( 6+rsub ) = zero
597 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
598 $ result( 6+rsub ) = ulpinv
603 IF( isort.EQ.1 )
THEN
607 IF(
sslect( wr( i ), wi( i ) ) .OR.
608 $
sslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
610 IF( (
sslect( wr( i+1 ), wi( i+1 ) ) .OR.
611 $
sslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
612 $ ( .NOT.(
sslect( wr( i ),
613 $ wi( i ) ) .OR.
sslect( wr( i ),
614 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
619 $ result( 13 ) = ulpinv
627 IF( lwork.GE.n+( n*n ) / 2 )
THEN
634 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
635 CALL
sgeesx(
'V', sort,
sslect,
'B', n, ht, lda, sdim1, wrt,
636 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
637 $ iwork, liwork, bwork, iinfo )
638 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
639 result( 14 ) = ulpinv
640 result( 15 ) = ulpinv
641 IF( jtype.NE.22 )
THEN
642 WRITE( nounit, fmt = 9998 )
'SGEESX3', iinfo, n, jtype,
645 WRITE( nounit, fmt = 9999 )
'SGEESX3', iinfo, n,
655 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
656 $ result( 10 ) = ulpinv
658 IF( h( i, j ).NE.ht( i, j ) )
659 $ result( 11 ) = ulpinv
660 IF( vs( i, j ).NE.vs1( i, j ) )
661 $ result( 12 ) = ulpinv
665 $ result( 13 ) = ulpinv
669 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
670 CALL
sgeesx(
'N', sort,
sslect,
'B', n, ht, lda, sdim1, wrt,
671 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
672 $ iwork, liwork, bwork, iinfo )
673 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
674 result( 14 ) = ulpinv
675 result( 15 ) = ulpinv
676 IF( jtype.NE.22 )
THEN
677 WRITE( nounit, fmt = 9998 )
'SGEESX4', iinfo, n, jtype,
680 WRITE( nounit, fmt = 9999 )
'SGEESX4', iinfo, n,
689 IF( rcnde1.NE.rconde )
690 $ result( 14 ) = ulpinv
691 IF( rcndv1.NE.rcondv )
692 $ result( 15 ) = ulpinv
697 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
698 $ result( 10 ) = ulpinv
700 IF( h( i, j ).NE.ht( i, j ) )
701 $ result( 11 ) = ulpinv
702 IF( vs( i, j ).NE.vs1( i, j ) )
703 $ result( 12 ) = ulpinv
707 $ result( 13 ) = ulpinv
711 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
712 CALL
sgeesx(
'V', sort,
sslect,
'E', n, ht, lda, sdim1, wrt,
713 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
714 $ iwork, liwork, bwork, iinfo )
715 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
716 result( 14 ) = ulpinv
717 IF( jtype.NE.22 )
THEN
718 WRITE( nounit, fmt = 9998 )
'SGEESX5', iinfo, n, jtype,
721 WRITE( nounit, fmt = 9999 )
'SGEESX5', iinfo, n,
730 IF( rcnde1.NE.rconde )
731 $ result( 14 ) = ulpinv
736 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
737 $ result( 10 ) = ulpinv
739 IF( h( i, j ).NE.ht( i, j ) )
740 $ result( 11 ) = ulpinv
741 IF( vs( i, j ).NE.vs1( i, j ) )
742 $ result( 12 ) = ulpinv
746 $ result( 13 ) = ulpinv
750 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
751 CALL
sgeesx(
'N', sort,
sslect,
'E', n, ht, lda, sdim1, wrt,
752 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
753 $ iwork, liwork, bwork, iinfo )
754 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
755 result( 14 ) = ulpinv
756 IF( jtype.NE.22 )
THEN
757 WRITE( nounit, fmt = 9998 )
'SGEESX6', iinfo, n, jtype,
760 WRITE( nounit, fmt = 9999 )
'SGEESX6', iinfo, n,
769 IF( rcnde1.NE.rconde )
770 $ result( 14 ) = ulpinv
775 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
776 $ result( 10 ) = ulpinv
778 IF( h( i, j ).NE.ht( i, j ) )
779 $ result( 11 ) = ulpinv
780 IF( vs( i, j ).NE.vs1( i, j ) )
781 $ result( 12 ) = ulpinv
785 $ result( 13 ) = ulpinv
789 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
790 CALL
sgeesx(
'V', sort,
sslect,
'V', n, ht, lda, sdim1, wrt,
791 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
792 $ iwork, liwork, bwork, iinfo )
793 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
794 result( 15 ) = ulpinv
795 IF( jtype.NE.22 )
THEN
796 WRITE( nounit, fmt = 9998 )
'SGEESX7', iinfo, n, jtype,
799 WRITE( nounit, fmt = 9999 )
'SGEESX7', iinfo, n,
808 IF( rcndv1.NE.rcondv )
809 $ result( 15 ) = ulpinv
814 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
815 $ result( 10 ) = ulpinv
817 IF( h( i, j ).NE.ht( i, j ) )
818 $ result( 11 ) = ulpinv
819 IF( vs( i, j ).NE.vs1( i, j ) )
820 $ result( 12 ) = ulpinv
824 $ result( 13 ) = ulpinv
828 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
829 CALL
sgeesx(
'N', sort,
sslect,
'V', n, ht, lda, sdim1, wrt,
830 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
831 $ iwork, liwork, bwork, iinfo )
832 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
833 result( 15 ) = ulpinv
834 IF( jtype.NE.22 )
THEN
835 WRITE( nounit, fmt = 9998 )
'SGEESX8', iinfo, n, jtype,
838 WRITE( nounit, fmt = 9999 )
'SGEESX8', iinfo, n,
847 IF( rcndv1.NE.rcondv )
848 $ result( 15 ) = ulpinv
853 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
854 $ result( 10 ) = ulpinv
856 IF( h( i, j ).NE.ht( i, j ) )
857 $ result( 11 ) = ulpinv
858 IF( vs( i, j ).NE.vs1( i, j ) )
859 $ result( 12 ) = ulpinv
863 $ result( 13 ) = ulpinv
880 eps = max( ulp, epsin )
883 selval( i ) = .false.
884 selwr( i ) = wrtmp( i )
885 selwi( i ) = witmp( i )
892 IF( wrtmp( j ).LT.vrmin )
THEN
898 wrtmp( kmin ) = wrtmp( i )
899 witmp( kmin ) = witmp( i )
903 ipnt( i ) = ipnt( kmin )
907 selval( ipnt( islct( i ) ) ) = .true.
912 CALL
slacpy(
'F', n, n, a, lda, ht, lda )
913 CALL
sgeesx(
'N',
'S',
sslect,
'B', n, ht, lda, sdim1, wrt,
914 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
915 $ iwork, liwork, bwork, iinfo )
916 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
917 result( 16 ) = ulpinv
918 result( 17 ) = ulpinv
919 WRITE( nounit, fmt = 9999 )
'SGEESX9', iinfo, n, iseed( 1 )
927 anorm =
slange(
'1', n, n, a, lda, work )
928 v = max(
REAL( n )*eps*anorm, smlnum )
931 IF( v.GT.rcondv )
THEN
936 IF( v.GT.rcdvin )
THEN
941 tol = max( tol, smlnum / eps )
942 tolin = max( tolin, smlnum / eps )
943 IF( eps*( rcdein-tolin ).GT.rconde+tol )
THEN
944 result( 16 ) = ulpinv
945 ELSE IF( rcdein-tolin.GT.rconde+tol )
THEN
946 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
947 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) )
THEN
948 result( 16 ) = ulpinv
949 ELSE IF( rcdein+tolin.LT.rconde-tol )
THEN
950 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
958 IF( v.GT.rcondv*rconde )
THEN
963 IF( v.GT.rcdvin*rcdein )
THEN
968 tol = max( tol, smlnum / eps )
969 tolin = max( tolin, smlnum / eps )
970 IF( eps*( rcdvin-tolin ).GT.rcondv+tol )
THEN
971 result( 17 ) = ulpinv
972 ELSE IF( rcdvin-tolin.GT.rcondv+tol )
THEN
973 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
974 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) )
THEN
975 result( 17 ) = ulpinv
976 ELSE IF( rcdvin+tolin.LT.rcondv-tol )
THEN
977 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
986 9999 format(
' SGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
987 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
988 9998 format(
' SGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
989 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )