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,
'.', / 9
x,
'N=',
987 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
988 9998
FORMAT(
' SGET24: ', a,
' returned INFO=', i6,
'.', / 9
x,
'N=',
989 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
LOGICAL function sslect(ZR, ZI)
SSLECT
REAL function slamch(CMACH)
SLAMCH
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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