339 SUBROUTINE sget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
340 $ H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
341 $ LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
342 $ RESULT, WORK, LWORK, IWORK, BWORK, INFO )
350 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
351 REAL RCDEIN, RCDVIN, THRESH
355 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
356 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
357 $ result( 17 ), vs( ldvs, * ), vs1( ldvs, * ),
358 $ wi( * ), wit( * ), witmp( * ), work( * ),
359 $ wr( * ), wrt( * ), wrtmp( * )
366 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
368 parameter( epsin = 5.9605e-8 )
372 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
374 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
375 $ smlnum, tmp, tol, tolin, ulp, ulpinv, v, vimin,
383 REAL SELWI( 20 ), SELWR( 20 )
386 INTEGER SELDIM, SELOPT
389 COMMON / sslct / selopt, seldim, selval, selwr, selwi
394 EXTERNAL SSLECT, SLAMCH, SLANGE
400 INTRINSIC abs, max, min, real, sign, sqrt
407 IF( thresh.LT.zero )
THEN
409 ELSE IF( nounit.LE.0 )
THEN
411 ELSE IF( n.LT.0 )
THEN
413 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
415 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n )
THEN
417 ELSE IF( lwork.LT.3*n )
THEN
422 CALL xerbla(
'SGET24', -info )
437 smlnum = slamch(
'Safe minimum' )
438 ulp = slamch(
'Precision' )
446 IF( isort.EQ.0 )
THEN
456 CALL slacpy(
'F', n, n, a, lda, h, lda )
457 CALL sgeesx(
'V', sort, sslect,
'N', n, h, lda, sdim, wr, wi,
458 $ vs, ldvs, rconde, rcondv, work, lwork, iwork,
459 $ liwork, bwork, iinfo )
460 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
461 result( 1+rsub ) = ulpinv
462 IF( jtype.NE.22 )
THEN
463 WRITE( nounit, fmt = 9998 )
'SGEESX1', iinfo, n, jtype,
466 WRITE( nounit, fmt = 9999 )
'SGEESX1', iinfo, n,
472 IF( isort.EQ.0 )
THEN
473 CALL scopy( n, wr, 1, wrtmp, 1 )
474 CALL scopy( n, wi, 1, witmp, 1 )
479 result( 1+rsub ) = zero
482 IF( h( i, j ).NE.zero )
483 $ result( 1+rsub ) = ulpinv
487 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
488 $ result( 1+rsub ) = ulpinv
491 IF( h( i+1, i ).NE.zero )
THEN
492 IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
493 $ zero .OR. sign( one, h( i+1, i ) ).EQ.
494 $ sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
502 CALL slacpy(
' ', n, n, a, lda, vs1, ldvs )
506 CALL sgemm(
'No transpose',
'No transpose', n, n, n, one, vs,
507 $ ldvs, h, lda, zero, ht, lda )
511 CALL sgemm(
'No transpose',
'Transpose', n, n, n, -one, ht,
512 $ lda, vs, ldvs, one, vs1, ldvs )
514 anorm = max( slange(
'1', n, n, a, lda, work ), smlnum )
515 wnorm = slange(
'1', n, n, vs1, ldvs, work )
517 IF( anorm.GT.wnorm )
THEN
518 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
520 IF( anorm.LT.one )
THEN
521 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
524 result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
531 CALL sort01(
'Columns', n, n, vs, ldvs, work, lwork,
536 result( 4+rsub ) = zero
538 IF( h( i, i ).NE.wr( i ) )
539 $ result( 4+rsub ) = ulpinv
542 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
543 $ result( 4+rsub ) = ulpinv
544 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
545 $ result( 4+rsub ) = ulpinv
548 IF( h( i+1, i ).NE.zero )
THEN
549 tmp = sqrt( abs( h( i+1, i ) ) )*
550 $ sqrt( abs( h( i, i+1 ) ) )
551 result( 4+rsub ) = max( result( 4+rsub ),
552 $ abs( wi( i )-tmp ) /
553 $ max( ulp*tmp, smlnum ) )
554 result( 4+rsub ) = max( result( 4+rsub ),
555 $ abs( wi( i+1 )+tmp ) /
556 $ max( ulp*tmp, smlnum ) )
557 ELSE IF( i.GT.1 )
THEN
558 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
559 $ wi( i ).NE.zero )result( 4+rsub ) = ulpinv
565 CALL slacpy(
'F', n, n, a, lda, ht, lda )
566 CALL sgeesx(
'N', sort, sslect,
'N', n, ht, lda, sdim, wrt,
567 $ wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
568 $ liwork, bwork, iinfo )
569 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
570 result( 5+rsub ) = ulpinv
571 IF( jtype.NE.22 )
THEN
572 WRITE( nounit, fmt = 9998 )
'SGEESX2', iinfo, n, jtype,
575 WRITE( nounit, fmt = 9999 )
'SGEESX2', iinfo, n,
582 result( 5+rsub ) = zero
585 IF( h( i, j ).NE.ht( i, j ) )
586 $ result( 5+rsub ) = ulpinv
592 result( 6+rsub ) = zero
594 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
595 $ result( 6+rsub ) = ulpinv
600 IF( isort.EQ.1 )
THEN
604 IF( sslect( wr( i ), wi( i ) ) .OR.
605 $ sslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
607 IF( ( sslect( wr( i+1 ), wi( i+1 ) ) .OR.
608 $ sslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
609 $ ( .NOT.( sslect( wr( i ),
610 $ wi( i ) ) .OR. sslect( wr( i ),
611 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
616 $ result( 13 ) = ulpinv
624 IF( lwork.GE.n+( n*n ) / 2 )
THEN
631 CALL slacpy(
'F', n, n, a, lda, ht, lda )
632 CALL sgeesx(
'V', sort, sslect,
'B', n, ht, lda, sdim1, wrt,
633 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
634 $ iwork, liwork, bwork, iinfo )
635 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
636 result( 14 ) = ulpinv
637 result( 15 ) = ulpinv
638 IF( jtype.NE.22 )
THEN
639 WRITE( nounit, fmt = 9998 )
'SGEESX3', iinfo, n, jtype,
642 WRITE( nounit, fmt = 9999 )
'SGEESX3', iinfo, n,
652 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
653 $ result( 10 ) = ulpinv
655 IF( h( i, j ).NE.ht( i, j ) )
656 $ result( 11 ) = ulpinv
657 IF( vs( i, j ).NE.vs1( i, j ) )
658 $ result( 12 ) = ulpinv
662 $ result( 13 ) = ulpinv
666 CALL slacpy(
'F', n, n, a, lda, ht, lda )
667 CALL sgeesx(
'N', sort, sslect,
'B', n, ht, lda, sdim1, wrt,
668 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
669 $ iwork, liwork, bwork, iinfo )
670 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
671 result( 14 ) = ulpinv
672 result( 15 ) = ulpinv
673 IF( jtype.NE.22 )
THEN
674 WRITE( nounit, fmt = 9998 )
'SGEESX4', iinfo, n, jtype,
677 WRITE( nounit, fmt = 9999 )
'SGEESX4', iinfo, n,
686 IF( rcnde1.NE.rconde )
687 $ result( 14 ) = ulpinv
688 IF( rcndv1.NE.rcondv )
689 $ result( 15 ) = ulpinv
694 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
695 $ result( 10 ) = ulpinv
697 IF( h( i, j ).NE.ht( i, j ) )
698 $ result( 11 ) = ulpinv
699 IF( vs( i, j ).NE.vs1( i, j ) )
700 $ result( 12 ) = ulpinv
704 $ result( 13 ) = ulpinv
708 CALL slacpy(
'F', n, n, a, lda, ht, lda )
709 CALL sgeesx(
'V', sort, sslect,
'E', n, ht, lda, sdim1, wrt,
710 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
711 $ iwork, liwork, bwork, iinfo )
712 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
713 result( 14 ) = ulpinv
714 IF( jtype.NE.22 )
THEN
715 WRITE( nounit, fmt = 9998 )
'SGEESX5', iinfo, n, jtype,
718 WRITE( nounit, fmt = 9999 )
'SGEESX5', iinfo, n,
727 IF( rcnde1.NE.rconde )
728 $ result( 14 ) = ulpinv
733 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
734 $ result( 10 ) = ulpinv
736 IF( h( i, j ).NE.ht( i, j ) )
737 $ result( 11 ) = ulpinv
738 IF( vs( i, j ).NE.vs1( i, j ) )
739 $ result( 12 ) = ulpinv
743 $ result( 13 ) = ulpinv
747 CALL slacpy(
'F', n, n, a, lda, ht, lda )
748 CALL sgeesx(
'N', sort, sslect,
'E', n, ht, lda, sdim1, wrt,
749 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
750 $ iwork, liwork, bwork, iinfo )
751 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
752 result( 14 ) = ulpinv
753 IF( jtype.NE.22 )
THEN
754 WRITE( nounit, fmt = 9998 )
'SGEESX6', iinfo, n, jtype,
757 WRITE( nounit, fmt = 9999 )
'SGEESX6', iinfo, n,
766 IF( rcnde1.NE.rconde )
767 $ result( 14 ) = ulpinv
772 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
773 $ result( 10 ) = ulpinv
775 IF( h( i, j ).NE.ht( i, j ) )
776 $ result( 11 ) = ulpinv
777 IF( vs( i, j ).NE.vs1( i, j ) )
778 $ result( 12 ) = ulpinv
782 $ result( 13 ) = ulpinv
786 CALL slacpy(
'F', n, n, a, lda, ht, lda )
787 CALL sgeesx(
'V', sort, sslect,
'V', n, ht, lda, sdim1, wrt,
788 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
789 $ iwork, liwork, bwork, iinfo )
790 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
791 result( 15 ) = ulpinv
792 IF( jtype.NE.22 )
THEN
793 WRITE( nounit, fmt = 9998 )
'SGEESX7', iinfo, n, jtype,
796 WRITE( nounit, fmt = 9999 )
'SGEESX7', iinfo, n,
805 IF( rcndv1.NE.rcondv )
806 $ result( 15 ) = ulpinv
811 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
812 $ result( 10 ) = ulpinv
814 IF( h( i, j ).NE.ht( i, j ) )
815 $ result( 11 ) = ulpinv
816 IF( vs( i, j ).NE.vs1( i, j ) )
817 $ result( 12 ) = ulpinv
821 $ result( 13 ) = ulpinv
825 CALL slacpy(
'F', n, n, a, lda, ht, lda )
826 CALL sgeesx(
'N', sort, sslect,
'V', n, ht, lda, sdim1, wrt,
827 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
828 $ iwork, liwork, bwork, iinfo )
829 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
830 result( 15 ) = ulpinv
831 IF( jtype.NE.22 )
THEN
832 WRITE( nounit, fmt = 9998 )
'SGEESX8', iinfo, n, jtype,
835 WRITE( nounit, fmt = 9999 )
'SGEESX8', iinfo, n,
844 IF( rcndv1.NE.rcondv )
845 $ result( 15 ) = ulpinv
850 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
851 $ result( 10 ) = ulpinv
853 IF( h( i, j ).NE.ht( i, j ) )
854 $ result( 11 ) = ulpinv
855 IF( vs( i, j ).NE.vs1( i, j ) )
856 $ result( 12 ) = ulpinv
860 $ result( 13 ) = ulpinv
877 eps = max( ulp, epsin )
880 selval( i ) = .false.
881 selwr( i ) = wrtmp( i )
882 selwi( i ) = witmp( i )
889 IF( wrtmp( j ).LT.vrmin )
THEN
895 wrtmp( kmin ) = wrtmp( i )
896 witmp( kmin ) = witmp( i )
900 ipnt( i ) = ipnt( kmin )
904 selval( ipnt( islct( i ) ) ) = .true.
909 CALL slacpy(
'F', n, n, a, lda, ht, lda )
910 CALL sgeesx(
'N',
'S', sslect,
'B', n, ht, lda, sdim1, wrt,
911 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
912 $ iwork, liwork, bwork, iinfo )
913 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
914 result( 16 ) = ulpinv
915 result( 17 ) = ulpinv
916 WRITE( nounit, fmt = 9999 )
'SGEESX9', iinfo, n, iseed( 1 )
924 anorm = slange(
'1', n, n, a, lda, work )
925 v = max( real( n )*eps*anorm, smlnum )
928 IF( v.GT.rcondv )
THEN
933 IF( v.GT.rcdvin )
THEN
938 tol = max( tol, smlnum / eps )
939 tolin = max( tolin, smlnum / eps )
940 IF( eps*( rcdein-tolin ).GT.rconde+tol )
THEN
941 result( 16 ) = ulpinv
942 ELSE IF( rcdein-tolin.GT.rconde+tol )
THEN
943 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
944 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) )
THEN
945 result( 16 ) = ulpinv
946 ELSE IF( rcdein+tolin.LT.rconde-tol )
THEN
947 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
955 IF( v.GT.rcondv*rconde )
THEN
960 IF( v.GT.rcdvin*rcdein )
THEN
965 tol = max( tol, smlnum / eps )
966 tolin = max( tolin, smlnum / eps )
967 IF( eps*( rcdvin-tolin ).GT.rcondv+tol )
THEN
968 result( 17 ) = ulpinv
969 ELSE IF( rcdvin-tolin.GT.rcondv+tol )
THEN
970 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
971 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) )
THEN
972 result( 17 ) = ulpinv
973 ELSE IF( rcdvin+tolin.LT.rcondv-tol )
THEN
974 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
983 9999
FORMAT(
' SGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
984 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
985 9998
FORMAT(
' SGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
986 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
subroutine xerbla(srname, info)
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 sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
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 sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01