sget38.f

Go to the documentation of this file.

*> \brief \b SGET38
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SGET38( RMAX, LMAX, NINFO, KNT, NIN )
*
*       .. Scalar Arguments ..
*       INTEGER            KNT, NIN
*       ..
*       .. Array Arguments ..
*       INTEGER            LMAX( 3 ), NINFO( 3 )
*       REAL               RMAX( 3 )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SGET38 tests STRSEN, a routine for estimating condition numbers of a
*> cluster of eigenvalues and/or its associated right invariant subspace
*>
*> The test matrices are read from a file with logical unit number NIN.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[out] RMAX
*> \verbatim
*>          RMAX is REAL array, dimension (3)
*>          Values of the largest test ratios.
*>          RMAX(1) = largest residuals from SHST01 or comparing
*>                    different calls to STRSEN
*>          RMAX(2) = largest error in reciprocal condition
*>                    numbers taking their conditioning into account
*>          RMAX(3) = largest error in reciprocal condition
*>                    numbers not taking their conditioning into
*>                    account (may be larger than RMAX(2))
*> \endverbatim
*>
*> \param[out] LMAX
*> \verbatim
*>          LMAX is INTEGER array, dimension (3)
*>          LMAX(i) is example number where largest test ratio
*>          RMAX(i) is achieved. Also:
*>          If SGEHRD returns INFO nonzero on example i, LMAX(1)=i
*>          If SHSEQR returns INFO nonzero on example i, LMAX(2)=i
*>          If STRSEN returns INFO nonzero on example i, LMAX(3)=i
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER array, dimension (3)
*>          NINFO(1) = No. of times SGEHRD returned INFO nonzero
*>          NINFO(2) = No. of times SHSEQR returned INFO nonzero
*>          NINFO(3) = No. of times STRSEN returned INFO nonzero
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*>          KNT is INTEGER
*>          Total number of examples tested.
*> \endverbatim
*>
*> \param[in] NIN
*> \verbatim
*>          NIN is INTEGER
*>          Input logical unit number.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup single_eig
*
*  =====================================================================
      SUBROUTINE sget38( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            KNT, NIN
*     ..
*     .. Array Arguments ..
      INTEGER            LMAX( 3 ), NINFO( 3 )
      REAL               RMAX( 3 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
      REAL               EPSIN
      parameter( epsin = 5.9605e-8 )
      INTEGER            LDT, LWORK
      parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
      INTEGER            LIWORK
      parameter( liwork = ldt*ldt )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, ISCL, ITMP, J, KMIN, M, N, NDIM
      REAL               BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
     $                   SMLNUM, STMP, TNRM, TOL, TOLIN, V, VIMIN, VMAX,
     $                   VMUL, VRMIN
*     ..
*     .. Local Arrays ..
      LOGICAL            SELECT( LDT )
      INTEGER            IPNT( LDT ), ISELEC( LDT ), IWORK( LIWORK )
      REAL               Q( LDT, LDT ), QSAV( LDT, LDT ),
     $                   QTMP( LDT, LDT ), RESULT( 2 ), T( LDT, LDT ),
     $                   TMP( LDT, LDT ), TSAV( LDT, LDT ),
     $                   TSAV1( LDT, LDT ), TTMP( LDT, LDT ), VAL( 3 ),
     $                   WI( LDT ), WITMP( LDT ), WORK( LWORK ),
     $                   WR( LDT ), WRTMP( LDT )
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLANGE
      EXTERNAL           slamch, slange
*     ..
*     .. External Subroutines ..
      EXTERNAL           scopy, sgehrd, shseqr, shst01, slabad, slacpy,
     $                   sorghr, sscal, strsen
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, real, sqrt
*     ..
*     .. Executable Statements ..
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' ) / eps
      bignum = one / smlnum
      CALL slabad( smlnum, bignum )
*
*     EPSIN = 2**(-24) = precision to which input data computed
*
      eps = max( eps, epsin )
      rmax( 1 ) = zero
      rmax( 2 ) = zero
      rmax( 3 ) = zero
      lmax( 1 ) = 0
      lmax( 2 ) = 0
      lmax( 3 ) = 0
      knt = 0
      ninfo( 1 ) = 0
      ninfo( 2 ) = 0
      ninfo( 3 ) = 0
*
      val( 1 ) = sqrt( smlnum )
      val( 2 ) = one
      val( 3 ) = sqrt( sqrt( bignum ) )
*
*     Read input data until N=0.  Assume input eigenvalues are sorted
*     lexicographically (increasing by real part, then decreasing by
*     imaginary part)
*
   10 CONTINUE
      READ( nin, fmt = * )n, ndim
      IF( n.EQ.0 )
     $   RETURN
      READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
      DO 20 i = 1, n
         READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
   20 CONTINUE
      READ( nin, fmt = * )sin, sepin
*
      tnrm = slange( 'M', n, n, tmp, ldt, work )
      DO 160 iscl = 1, 3
*
*        Scale input matrix
*
         knt = knt + 1
         CALL slacpy( 'F', n, n, tmp, ldt, t, ldt )
         vmul = val( iscl )
         DO 30 i = 1, n
            CALL sscal( n, vmul, t( 1, i ), 1 )
   30    CONTINUE
         IF( tnrm.EQ.zero )
     $      vmul = one
         CALL slacpy( 'F', n, n, t, ldt, tsav, ldt )
*
*        Compute Schur form
*
         CALL sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
     $                info )
         IF( info.NE.0 ) THEN
            lmax( 1 ) = knt
            ninfo( 1 ) = ninfo( 1 ) + 1
            GO TO 160
         END IF
*
*        Generate orthogonal matrix
*
         CALL slacpy( 'L', n, n, t, ldt, q, ldt )
         CALL sorghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
     $                info )
*
*        Compute Schur form
*
         CALL shseqr( 'S', 'V', n, 1, n, t, ldt, wr, wi, q, ldt, work,
     $                lwork, info )
         IF( info.NE.0 ) THEN
            lmax( 2 ) = knt
            ninfo( 2 ) = ninfo( 2 ) + 1
            GO TO 160
         END IF
*
*        Sort, select eigenvalues
*
         DO 40 i = 1, n
            ipnt( i ) = i
            SELECT( i ) = .false.
   40    CONTINUE
         CALL scopy( n, wr, 1, wrtmp, 1 )
         CALL scopy( n, wi, 1, witmp, 1 )
         DO 60 i = 1, n - 1
            kmin = i
            vrmin = wrtmp( i )
            vimin = witmp( i )
            DO 50 j = i + 1, n
               IF( wrtmp( j ).LT.vrmin ) THEN
                  kmin = j
                  vrmin = wrtmp( j )
                  vimin = witmp( j )
               END IF
   50       CONTINUE
            wrtmp( kmin ) = wrtmp( i )
            witmp( kmin ) = witmp( i )
            wrtmp( i ) = vrmin
            witmp( i ) = vimin
            itmp = ipnt( i )
            ipnt( i ) = ipnt( kmin )
            ipnt( kmin ) = itmp
   60    CONTINUE
         DO 70 i = 1, ndim
            SELECT( ipnt( iselec( i ) ) ) = .true.
   70    CONTINUE
*
*        Compute condition numbers
*
         CALL slacpy( 'F', n, n, q, ldt, qsav, ldt )
         CALL slacpy( 'F', n, n, t, ldt, tsav1, ldt )
         CALL strsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wrtmp, witmp,
     $                m, s, sep, work, lwork, iwork, liwork, info )
         IF( info.NE.0 ) THEN
            lmax( 3 ) = knt
            ninfo( 3 ) = ninfo( 3 ) + 1
            GO TO 160
         END IF
         septmp = sep / vmul
         stmp = s
*
*        Compute residuals
*
         CALL shst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
     $                result )
         vmax = max( result( 1 ), result( 2 ) )
         IF( vmax.GT.rmax( 1 ) ) THEN
            rmax( 1 ) = vmax
            IF( ninfo( 1 ).EQ.0 )
     $         lmax( 1 ) = knt
         END IF
*
*        Compare condition number for eigenvalue cluster
*        taking its condition number into account
*
         v = max( two*real( n )*eps*tnrm, smlnum )
         IF( tnrm.EQ.zero )
     $      v = one
         IF( v.GT.septmp ) THEN
            tol = one
         ELSE
            tol = v / septmp
         END IF
         IF( v.GT.sepin ) THEN
            tolin = one
         ELSE
            tolin = v / sepin
         END IF
         tol = max( tol, smlnum / eps )
         tolin = max( tolin, smlnum / eps )
         IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
            vmax = one / eps
         ELSE IF( sin-tolin.GT.stmp+tol ) THEN
            vmax = ( sin-tolin ) / ( stmp+tol )
         ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
            vmax = one / eps
         ELSE IF( sin+tolin.LT.stmp-tol ) THEN
            vmax = ( stmp-tol ) / ( sin+tolin )
         ELSE
            vmax = one
         END IF
         IF( vmax.GT.rmax( 2 ) ) THEN
            rmax( 2 ) = vmax
            IF( ninfo( 2 ).EQ.0 )
     $         lmax( 2 ) = knt
         END IF
*
*        Compare condition numbers for invariant subspace
*        taking its condition number into account
*
         IF( v.GT.septmp*stmp ) THEN
            tol = septmp
         ELSE
            tol = v / stmp
         END IF
         IF( v.GT.sepin*sin ) THEN
            tolin = sepin
         ELSE
            tolin = v / sin
         END IF
         tol = max( tol, smlnum / eps )
         tolin = max( tolin, smlnum / eps )
         IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
            vmax = one / eps
         ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
            vmax = ( sepin-tolin ) / ( septmp+tol )
         ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
            vmax = one / eps
         ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
            vmax = ( septmp-tol ) / ( sepin+tolin )
         ELSE
            vmax = one
         END IF
         IF( vmax.GT.rmax( 2 ) ) THEN
            rmax( 2 ) = vmax
            IF( ninfo( 2 ).EQ.0 )
     $         lmax( 2 ) = knt
         END IF
*
*        Compare condition number for eigenvalue cluster
*        without taking its condition number into account
*
         IF( sin.LE.real( 2*n )*eps .AND. stmp.LE.real( 2*n )*eps ) THEN
            vmax = one
         ELSE IF( eps*sin.GT.stmp ) THEN
            vmax = one / eps
         ELSE IF( sin.GT.stmp ) THEN
            vmax = sin / stmp
         ELSE IF( sin.LT.eps*stmp ) THEN
            vmax = one / eps
         ELSE IF( sin.LT.stmp ) THEN
            vmax = stmp / sin
         ELSE
            vmax = one
         END IF
         IF( vmax.GT.rmax( 3 ) ) THEN
            rmax( 3 ) = vmax
            IF( ninfo( 3 ).EQ.0 )
     $         lmax( 3 ) = knt
         END IF
*
*        Compare condition numbers for invariant subspace
*        without taking its condition number into account
*
         IF( sepin.LE.v .AND. septmp.LE.v ) THEN
            vmax = one
         ELSE IF( eps*sepin.GT.septmp ) THEN
            vmax = one / eps
         ELSE IF( sepin.GT.septmp ) THEN
            vmax = sepin / septmp
         ELSE IF( sepin.LT.eps*septmp ) THEN
            vmax = one / eps
         ELSE IF( sepin.LT.septmp ) THEN
            vmax = septmp / sepin
         ELSE
            vmax = one
         END IF
         IF( vmax.GT.rmax( 3 ) ) THEN
            rmax( 3 ) = vmax
            IF( ninfo( 3 ).EQ.0 )
     $         lmax( 3 ) = knt
         END IF
*
*        Compute eigenvalue condition number only and compare
*        Update Q
*
         vmax = zero
         CALL slacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
         CALL slacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
         septmp = -one
         stmp = -one
         CALL strsen( 'e', 'v', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
     $                WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
     $                LIWORK, INFO )
.NE.         IF( INFO0 ) THEN
            LMAX( 3 ) = KNT
            NINFO( 3 ) = NINFO( 3 ) + 1
            GO TO 160
         END IF
.NE.         IF( SSTMP )
     $      VMAX = ONE / EPS
.NE.         IF( -ONESEPTMP )
     $      VMAX = ONE / EPS
         DO 90 I = 1, N
            DO 80 J = 1, N
.NE.               IF( TTMP( I, J )T( I, J ) )
     $            VMAX = ONE / EPS
.NE.               IF( QTMP( I, J )Q( I, J ) )
     $            VMAX = ONE / EPS
   80       CONTINUE
   90    CONTINUE
*
*        Compute invariant subspace condition number only and compare
*        Update Q
*
         CALL SLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
         CALL SLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
         SEPTMP = -ONE
         STMP = -ONE
         CALL STRSEN( 'v', 'v', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
     $                WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
     $                LIWORK, INFO )
.NE.         IF( INFO0 ) THEN
            LMAX( 3 ) = KNT
            NINFO( 3 ) = NINFO( 3 ) + 1
            GO TO 160
         END IF
.NE.         IF( -ONESTMP )
     $      VMAX = ONE / EPS
.NE.         IF( SEPSEPTMP )
     $      VMAX = ONE / EPS
         DO 110 I = 1, N
            DO 100 J = 1, N
.NE.               IF( TTMP( I, J )T( I, J ) )
     $            VMAX = ONE / EPS
.NE.               IF( QTMP( I, J )Q( I, J ) )
     $            VMAX = ONE / EPS
  100       CONTINUE
  110    CONTINUE
*
*        Compute eigenvalue condition number only and compare
*        Do not update Q
*
         CALL SLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
         CALL SLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
         SEPTMP = -ONE
         STMP = -ONE
         CALL STRSEN( 'e', 'n', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
     $                WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
     $                LIWORK, INFO )
.NE.         IF( INFO0 ) THEN
            LMAX( 3 ) = KNT
            NINFO( 3 ) = NINFO( 3 ) + 1
            GO TO 160
         END IF
.NE.         IF( SSTMP )
     $      VMAX = ONE / EPS
.NE.         IF( -ONESEPTMP )
     $      VMAX = ONE / EPS
         DO 130 I = 1, N
            DO 120 J = 1, N
.NE.               IF( TTMP( I, J )T( I, J ) )
     $            VMAX = ONE / EPS
.NE.               IF( QTMP( I, J )QSAV( I, J ) )
     $            VMAX = ONE / EPS
  120       CONTINUE
  130    CONTINUE
*
*        Compute invariant subspace condition number only and compare
*        Do not update Q
*
         CALL SLACPY( 'f', N, N, TSAV1, LDT, TTMP, LDT )
         CALL SLACPY( 'f', N, N, QSAV, LDT, QTMP, LDT )
         SEPTMP = -ONE
         STMP = -ONE
         CALL STRSEN( 'v', 'n', SELECT, N, TTMP, LDT, QTMP, LDT, WRTMP,
     $                WITMP, M, STMP, SEPTMP, WORK, LWORK, IWORK,
     $                LIWORK, INFO )
.NE.         IF( INFO0 ) THEN
            LMAX( 3 ) = KNT
            NINFO( 3 ) = NINFO( 3 ) + 1
            GO TO 160
         END IF
.NE.         IF( -ONESTMP )
     $      VMAX = ONE / EPS
.NE.         IF( SEPSEPTMP )
     $      VMAX = ONE / EPS
         DO 150 I = 1, N
            DO 140 J = 1, N
.NE.               IF( TTMP( I, J )T( I, J ) )
     $            VMAX = ONE / EPS
.NE.               IF( QTMP( I, J )QSAV( I, J ) )
     $            VMAX = ONE / EPS
  140       CONTINUE
  150    CONTINUE
.GT.         IF( VMAXRMAX( 1 ) ) THEN
            RMAX( 1 ) = VMAX
.EQ.            IF( NINFO( 1 )0 )
     $         LMAX( 1 ) = KNT
         END IF
  160 CONTINUE
      GO TO 10
*
*     End of SGET38
*
      END