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
397 EXTERNAL sslect, slamch, slange
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,
')' )
subroutine sgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, IWORK, BWORK, INFO)
SGET24
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY