dget40.f

Go to the documentation of this file.

*> \brief \b DGET40
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET40( RMAX, LMAX, NINFO, KNT, NIN )
*
*       .. Scalar Arguments ..
*      INTEGER            KNT, LMAX, NIN
*      DOUBLE PRECISION   RMAX
*      ..
*       .. Array Arguments ..
*      INTEGER            NINFO( 2 )
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGET40 tests DTGEXC, a routine for swapping adjacent blocks (either
*> 1 by 1 or 2 by 2) on the diagonal of a pencil in real generalized Schur form.
*> Thus, DTGEXC computes an orthogonal matrices Q and Z such that
*>
*>     Q' * ( [ A B ], [ D E ] ) * Z  = ( [ C1 B1 ], [ F1 E1 ] )
*>          ( [ 0 C ]  [   F ] )        ( [ 0  A1 ]  [    D1]  )
*>
*> where (C1,F1) is similar to (C,F) and (A1,D1) is similar to (A,D).
*> Both (A,D) and (C,F) are assumed to be in standard form
*> and (A1,D1) and (C1,F1) are returned with the
*> same properties.
*> \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 (2)
*>          NINFO( 1 ) = DTGEXC without accumulation returned INFO nonzero
*>          NINFO( 2 ) = DTGEXC with accumulation 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 double_eig
*
*  =====================================================================
      SUBROUTINE dget40( 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, LMAX, NIN
      DOUBLE PRECISION   RMAX
*     ..
*     .. Array Arguments ..
      INTEGER            NINFO( 2 )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
      INTEGER            LDT, LWORK
      parameter( ldt = 10, lwork = 100 + 4*ldt + 16 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1,
     $                   ILST2, ILSTSV, J, LOC, N
      DOUBLE PRECISION   EPS, RES
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   Q( LDT, LDT ), Z( LDT, LDT ), RESULT( 4 ),
     $                   T( LDT, LDT ), T1( LDT, LDT ), T2( LDT, LDT ),
     $                   S( LDT, LDT ), S1( LDT, LDT ), S2( LDT, LDT ),
     $                   TMP( LDT, LDT ), WORK( LWORK )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dhst01, dlacpy, dlaset, dtgexc
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, sign
*     ..
*     .. Executable Statements ..
*
      eps = dlamch( 'P' )
      rmax = zero
      lmax = 0
      knt = 0
      ninfo( 1 ) = 0
      ninfo( 2 ) = 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, t, ldt )
      CALL dlacpy( 'F', n, n, tmp, ldt, t1, ldt )
      CALL dlacpy( 'F', n, n, tmp, ldt, t2, ldt )
      DO 25 i = 1, n
         READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
   25 CONTINUE
      CALL dlacpy( 'F', n, n, tmp, ldt, s, ldt )
      CALL dlacpy( 'F', n, n, tmp, ldt, s1, ldt )
      CALL dlacpy( 'F', n, n, tmp, ldt, s2, ldt )
      ifstsv = ifst
      ilstsv = ilst
      ifst1 = ifst
      ilst1 = ilst
      ifst2 = ifst
      ilst2 = ilst
      res = zero
*
*     Test without accumulating Q and Z
*
      CALL dlaset( 'Full', n, n, zero, one, q, ldt )
      CALL dlaset( 'Full', n, n, zero, one, z, ldt )
      CALL dtgexc( .false., .false., n, t1, ldt, s1, ldt, q, ldt,
     $             z, ldt, ifst1, ilst1, work, lwork, ninfo( 1 ) )
      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
            IF( i.EQ.j .AND. z( i, j ).NE.one )
     $         res = res + one / eps
            IF( i.NE.j .AND. z( 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 dlaset( 'Full', n, n, zero, one, z, ldt )
      CALL dtgexc( .true., .true., n, t2, ldt, s2, ldt, q, ldt,
     $             z, ldt, ifst2, ilst2, work, lwork, ninfo( 2 ) )
*
*     Compare T1 with T2 and S1 with S2
*
      DO 60 i = 1, n
         DO 50 j = 1, n
            IF( t1( i, j ).NE.t2( i, j ) )
     $         res = res + one / eps
            IF( s1( i, j ).NE.s2( 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( ninfo( 1 ).NE.ninfo( 2 ) )
     $   res = res + one / eps
*
*     Test orthogonality of Q and Z and backward error on T2 and S2
*
      CALL dget51( 1, n, t, ldt, t2, ldt, q, ldt, z, ldt, work,
     $             result( 1 ) )
      CALL dget51( 1, n, s, ldt, s2, ldt, q, ldt, z, ldt, work,
     $             result( 2 ) )
      CALL dget51( 3, n, t, ldt, t2, ldt, q, ldt, q, ldt, work,
     $             result( 3 ) )
      CALL dget51( 3, n, t, ldt, t2, ldt, z, ldt, z, ldt, work,
     $             result( 4 ) )
      res = res + result( 1 ) + result( 2 ) + result( 3 ) + result( 4 )
*
*     Read next matrix pair
*
      GO TO 10
*
*     End of DGET40
*
      END