dget36.f

Go to the documentation of this file.

*> \brief \b DGET36
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET36( RMAX, LMAX, NINFO, KNT, NIN )
* 
*       .. Scalar Arguments ..
*       INTEGER            KNT, LMAX, NIN
*       DOUBLE PRECISION   RMAX
*       ..
*       .. Array Arguments ..
*       INTEGER            NINFO( 3 )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGET36 tests DTREXC, a routine for moving blocks (either 1 by 1 or
*> 2 by 2) on the diagonal of a matrix in real Schur form.  Thus, DLAEXC
*> computes an orthogonal matrix Q such that
*>
*>    Q' * T1 * Q  = T2
*>
*> and where one of the diagonal blocks of T1 (the one at row IFST) has
*> been moved to position ILST.
*>
*> The test code verifies that the residual Q'*T1*Q-T2 is small, that T2
*> is in Schur form, and that the final position of the IFST block is
*> ILST (within +-1).
*>
*> The test matrices are read from a file with logical unit number NIN.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[out] RMAX
*> \verbatim
*>          RMAX is DOUBLE PRECISION
*>          Value of the largest test ratio.
*> \endverbatim
*>
*> \param[out] LMAX
*> \verbatim
*>          LMAX is INTEGER
*>          Example number where largest test ratio achieved.
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER array, dimension (3)
*>          NINFO(J) is the number of examples where INFO=J.
*> \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. 
*
*> \date November 2011
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dget36( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            knt, lmax, nin
      DOUBLE PRECISION   rmax
*     ..
*     .. Array Arguments ..
      INTEGER            ninfo( 3 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, one
      parameter( zero = 0.0d0, one = 1.0d0 )
      INTEGER            ldt, lwork
      parameter( ldt = 10, lwork = 2*ldt*ldt )
*     ..
*     .. Local Scalars ..
      INTEGER            i, ifst, ifst1, ifst2, ifstsv, ilst, ilst1,
     $                   ilst2, ilstsv, info1, info2, j, loc, n
      DOUBLE PRECISION   eps, res
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   q( ldt, ldt ), result( 2 ), t1( ldt, ldt ),
     $                   t2( ldt, ldt ), tmp( ldt, ldt ), work( lwork )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlamch
      EXTERNAL           dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dhst01, dlacpy, dlaset, dtrexc
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sign
*     ..
*     .. Executable Statements ..
*
      eps = dlamch( 'P' )
      rmax = zero
      lmax = 0
      knt = 0
      ninfo( 1 ) = 0
      ninfo( 2 ) = 0
      ninfo( 3 ) = 0
*
*     Read input data until N=0
*
   10 continue
      READ( nin, fmt = * )n, ifst, ilst
      IF( n.EQ.0 )
     $   return
      knt = knt + 1
      DO 20 i = 1, n
         READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
   20 continue
      CALL dlacpy( 'F', n, n, tmp, ldt, t1, ldt )
      CALL dlacpy( 'F', n, n, tmp, ldt, t2, ldt )
      ifstsv = ifst
      ilstsv = ilst
      ifst1 = ifst
      ilst1 = ilst
      ifst2 = ifst
      ilst2 = ilst
      res = zero
*
*     Test without accumulating Q
*
      CALL dlaset( 'Full', n, n, zero, one, q, ldt )
      CALL dtrexc( 'N', n, t1, ldt, q, ldt, ifst1, ilst1, work, info1 )
      DO 40 i = 1, n
         DO 30 j = 1, n
            IF( i.EQ.j .AND. q( i, j ).NE.one )
     $         res = res + one / eps
            IF( i.NE.j .AND. q( i, j ).NE.zero )
     $         res = res + one / eps
   30    continue
   40 continue
*
*     Test with accumulating Q
*
      CALL dlaset( 'Full', n, n, zero, one, q, ldt )
      CALL dtrexc( 'V', n, t2, ldt, q, ldt, ifst2, ilst2, work, info2 )
*
*     Compare T1 with T2
*
      DO 60 i = 1, n
         DO 50 j = 1, n
            IF( t1( i, j ).NE.t2( i, j ) )
     $         res = res + one / eps
   50    continue
   60 continue
      IF( ifst1.NE.ifst2 )
     $   res = res + one / eps
      IF( ilst1.NE.ilst2 )
     $   res = res + one / eps
      IF( info1.NE.info2 )
     $   res = res + one / eps
*
*     Test for successful reordering of T2
*
      IF( info2.NE.0 ) THEN
         ninfo( info2 ) = ninfo( info2 ) + 1
      ELSE
         IF( abs( ifst2-ifstsv ).GT.1 )
     $      res = res + one / eps
         IF( abs( ilst2-ilstsv ).GT.1 )
     $      res = res + one / eps
      END IF
*
*     Test for small residual, and orthogonality of Q
*
      CALL dhst01( n, 1, n, tmp, ldt, t2, ldt, q, ldt, work, lwork,
     $             result )
      res = res + result( 1 ) + result( 2 )
*
*     Test for T2 being in Schur form
*
      loc = 1
   70 continue
      IF( t2( loc+1, loc ).NE.zero ) THEN
*
*        2 by 2 block
*
         IF( t2( loc, loc+1 ).EQ.zero .OR. t2( loc, loc ).NE.
     $       t2( loc+1, loc+1 ) .OR. sign( one, t2( loc, loc+1 ) ).EQ.
     $       sign( one, t2( loc+1, loc ) ) )res = res + one / eps
         DO 80 i = loc + 2, n
            IF( t2( i, loc ).NE.zero )
     $         res = res + one / res
            IF( t2( i, loc+1 ).NE.zero )
     $         res = res + one / res
   80    continue
         loc = loc + 2
      ELSE
*
*        1 by 1 block
*
         DO 90 i = loc + 1, n
            IF( t2( i, loc ).NE.zero )
     $         res = res + one / res
   90    continue
         loc = loc + 1
      END IF
      IF( loc.LT.n )
     $   go to 70
      IF( res.GT.rmax ) THEN
         rmax = res
         lmax = knt
      END IF
      go to 10
*
*     End of DGET36
*
      END