text: *> \brief \b DGET39
  *
  *  =========== DOCUMENTATION ===========
  *
  * Online html documentation available at
  *            http://www.netlib.org/lapack/explore-html/
  *
  *  Definition:
  *  ===========
  *
  *       SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
  *
  *       .. Scalar Arguments ..
  *       INTEGER            KNT, LMAX, NINFO
  *       DOUBLE PRECISION   RMAX
  *       ..
  *
  *
  *> \par Purpose:
  *  =============
  *>
  *> \verbatim
  *>
  *> DGET39 tests DLAQTR, a routine for solving the real or
  *> special complex quasi upper triangular system
  *>
  *>      op(T)*p = scale*c,
  *> or
  *>      op(T + iB)*(p+iq) = scale*(c+id),
  *>
  *> in real arithmetic. T is upper quasi-triangular.
  *> If it is complex, then the first diagonal block of T must be
  *> 1 by 1, B has the special structure
  *>
  *>                B = [ b(1) b(2) ... b(n) ]
  *>                    [       w            ]
  *>                    [           w        ]
  *>                    [              .     ]
  *>                    [                 w  ]
  *>
  *> op(A) = A or A', where A' denotes the conjugate transpose of
  *> the matrix A.
  *>
  *> On input, X = [ c ].  On output, X = [ p ].
  *>               [ d ]                  [ q ]
  *>
  *> Scale is an output less than or equal to 1, chosen to avoid
  *> overflow in X.
  *> This subroutine is specially designed for the condition number
  *> estimation in the eigenproblem routine DTRSNA.
  *>
  *> The test code verifies that the following residual is order 1:
  *>
  *>      ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
  *>    -----------------------------------------
  *>        max(ulp*(||T||+||B||)*(||x1||+||x2||),
  *>            (||T||+||B||)*smlnum/ulp,
  *>            smlnum)
  *>
  *> (The (||T||+||B||)*smlnum/ulp term accounts for possible
  *>  (gradual or nongradual) underflow in x1 and x2.)
  *> \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
  *>          Number of examples where INFO is nonzero.
  *> \endverbatim
  *>
  *> \param[out] KNT
  *> \verbatim
  *>          KNT is INTEGER
  *>          Total number of examples tested.
  *> \endverbatim
  *
  *  Authors:
  *  ========
  *
  *> \author Univ. of Tennessee
  *> \author Univ. of California Berkeley
  *> \author Univ. of Colorado Denver
  *> \author NAG Ltd.
  *
  *> \date December 2016
  *
  *> \ingroup double_eig
  *
  *  =====================================================================
        SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
  *
  *  -- LAPACK test routine (version 3.7.0) --
  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  *     December 2016
  *
  *     .. Scalar Arguments ..
        INTEGER            KNT, LMAX, NINFO
        DOUBLE PRECISION   RMAX
  *     ..
  *
  *  =====================================================================
  *
  *     .. Parameters ..
        INTEGER            LDT, LDT2
        PARAMETER          ( LDT = 10, LDT2 = 2*LDT )
        DOUBLE PRECISION   ZERO, ONE
        PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  *     ..
  *     .. Local Scalars ..
        INTEGER            I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
       $                   NDIM
        DOUBLE PRECISION   BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
       $                   SCALE, SMLNUM, W, XNORM
  *     ..
  *     .. External Functions ..
        INTEGER            IDAMAX
        DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DLANGE
        EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
  *     ..
  *     .. External Subroutines ..
        EXTERNAL           DCOPY, DGEMV, DLABAD, DLAQTR
  *     ..
  *     .. Intrinsic Functions ..
        INTRINSIC          ABS, COS, DBLE, MAX, SIN, SQRT
  *     ..
  *     .. Local Arrays ..
        INTEGER            IDIM( 6 ), IVAL( 5, 5, 6 )
        DOUBLE PRECISION   B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
       $                   VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
       $                   VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
  *     ..
  *     .. Data statements ..
        DATA               IDIM / 4, 5*5 /
        DATA               IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
       $                   4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
       $                   0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
       $                   0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
       $                   4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
       $                   -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
       $                   5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
       $                   4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
  *     ..
  *     .. Executable Statements ..
  *
  *     Get machine parameters
  *
        EPS = DLAMCH( 'P' )
        SMLNUM = DLAMCH( 'S' )
        BIGNUM = ONE / SMLNUM
        CALL DLABAD( SMLNUM, BIGNUM )
  *
  *     Set up test case parameters
  *
        VM1( 1 ) = ONE
        VM1( 2 ) = SQRT( SMLNUM )
        VM1( 3 ) = SQRT( VM1( 2 ) )
        VM1( 4 ) = SQRT( BIGNUM )
        VM1( 5 ) = SQRT( VM1( 4 ) )
  *
        VM2( 1 ) = ONE
        VM2( 2 ) = SQRT( SMLNUM )
        VM2( 3 ) = SQRT( VM2( 2 ) )
        VM2( 4 ) = SQRT( BIGNUM )
        VM2( 5 ) = SQRT( VM2( 4 ) )
  *
        VM3( 1 ) = ONE
        VM3( 2 ) = SQRT( SMLNUM )
        VM3( 3 ) = SQRT( VM3( 2 ) )
        VM3( 4 ) = SQRT( BIGNUM )
        VM3( 5 ) = SQRT( VM3( 4 ) )
  *
        VM4( 1 ) = ONE
        VM4( 2 ) = SQRT( SMLNUM )
        VM4( 3 ) = SQRT( VM4( 2 ) )
        VM4( 4 ) = SQRT( BIGNUM )
        VM4( 5 ) = SQRT( VM4( 4 ) )
  *
        VM5( 1 ) = ONE
        VM5( 2 ) = EPS
        VM5( 3 ) = SQRT( SMLNUM )
  *
  *     Initialization
  *
        KNT = 0
        RMAX = ZERO
        NINFO = 0
        SMLNUM = SMLNUM / EPS
  *
  *     Begin test loop
  *
        DO 140 IVM5 = 1, 3
           DO 130 IVM4 = 1, 5
              DO 120 IVM3 = 1, 5
                 DO 110 IVM2 = 1, 5
                    DO 100 IVM1 = 1, 5
                       DO 90 NDIM = 1, 6
  *
                          N = IDIM( NDIM )
                          DO 20 I = 1, N
                             DO 10 J = 1, N
                                T( I, J ) = DBLE( IVAL( I, J, NDIM ) )*
       $                                    VM1( IVM1 )
                                IF( I.GE.J )
       $                           T( I, J ) = T( I, J )*VM5( IVM5 )
     10                      CONTINUE
     20                   CONTINUE
  *
                          W = ONE*VM2( IVM2 )
  *
                          DO 30 I = 1, N
                             B( I ) = COS( DBLE( I ) )*VM3( IVM3 )
     30                   CONTINUE
  *
                          DO 40 I = 1, 2*N
                             D( I ) = SIN( DBLE( I ) )*VM4( IVM4 )
     40                   CONTINUE
  *
                          NORM = DLANGE( '1', N, N, T, LDT, WORK )
                          K = IDAMAX( N, B, 1 )
                          NORMTB = NORM + ABS( B( K ) ) + ABS( W )
  *
                          CALL DCOPY( N, D, 1, X, 1 )
                          KNT = KNT + 1
                          CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
       $                               DUMM, SCALE, X, WORK, INFO )
                          IF( INFO.NE.0 )
       $                     NINFO = NINFO + 1
  *
  *                       || T*x - scale*d || /
  *                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
  *
                          CALL DCOPY( N, D, 1, Y, 1 )
                          CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
       $                              X, 1, -SCALE, Y, 1 )
                          XNORM = DASUM( N, X, 1 )
                          RESID = DASUM( N, Y, 1 )
                          DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
       $                          ( NORM*EPS )*XNORM )
                          RESID = RESID / DOMIN
                          IF( RESID.GT.RMAX ) THEN
                             RMAX = RESID
                             LMAX = KNT
                          END IF
  *
                          CALL DCOPY( N, D, 1, X, 1 )
                          KNT = KNT + 1
                          CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
       $                               DUMM, SCALE, X, WORK, INFO )
                          IF( INFO.NE.0 )
       $                     NINFO = NINFO + 1
  *
  *                       || T*x - scale*d || /
  *                         max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
  *
                          CALL DCOPY( N, D, 1, Y, 1 )
                          CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
       $                              1, -SCALE, Y, 1 )
                          XNORM = DASUM( N, X, 1 )
                          RESID = DASUM( N, Y, 1 )
                          DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
       $                          ( NORM*EPS )*XNORM )
                          RESID = RESID / DOMIN
                          IF( RESID.GT.RMAX ) THEN
                             RMAX = RESID
                             LMAX = KNT
                          END IF
  *
                          CALL DCOPY( 2*N, D, 1, X, 1 )
                          KNT = KNT + 1
                          CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
       $                               SCALE, X, WORK, INFO )
                          IF( INFO.NE.0 )
       $                     NINFO = NINFO + 1
  *
  *                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
  *                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
  *                                  smlnum/ulp * (||T||+||B||), smlnum )
  *
  *
                          CALL DCOPY( 2*N, D, 1, Y, 1 )
                          Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) +
       $                           SCALE*Y( 1 )
                          DO 50 I = 2, N
                             Y( I ) = W*X( I+N ) + SCALE*Y( I )
     50                   CONTINUE
                          CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
       $                              X, 1, -ONE, Y, 1 )
  *
                          Y( 1+N ) = DDOT( N, B, 1, X, 1 ) -
       $                             SCALE*Y( 1+N )
                          DO 60 I = 2, N
                             Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
     60                   CONTINUE
                          CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
       $                              X( 1+N ), 1, ONE, Y( 1+N ), 1 )
  *
                          RESID = DASUM( 2*N, Y, 1 )
                          DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
       $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
                          RESID = RESID / DOMIN
                          IF( RESID.GT.RMAX ) THEN
                             RMAX = RESID
                             LMAX = KNT
                          END IF
  *
                          CALL DCOPY( 2*N, D, 1, X, 1 )
                          KNT = KNT + 1
                          CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
       $                               SCALE, X, WORK, INFO )
                          IF( INFO.NE.0 )
       $                     NINFO = NINFO + 1
  *
  *                       ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
  *                          max(ulp*(||T||+||B||)*(||x1||+||x2||),
  *                                  smlnum/ulp * (||T||+||B||), smlnum )
  *
                          CALL DCOPY( 2*N, D, 1, Y, 1 )
                          Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
                          DO 70 I = 2, N
                             Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
       $                              SCALE*Y( I )
     70                   CONTINUE
                          CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
       $                              1, ONE, Y, 1 )
  *
                          Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
                          DO 80 I = 2, N
                             Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
       $                                SCALE*Y( I+N )
     80                   CONTINUE
                          CALL DGEMV( 'Transpose', N, N, ONE, T, LDT,
       $                              X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
  *
                          RESID = DASUM( 2*N, Y, 1 )
                          DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
       $                          EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
                          RESID = RESID / DOMIN
                          IF( RESID.GT.RMAX ) THEN
                             RMAX = RESID
                             LMAX = KNT
                          END IF
  *
     90                CONTINUE
    100             CONTINUE
    110          CONTINUE
    120       CONTINUE
    130    CONTINUE
    140 CONTINUE
  *
        RETURN
  *
  *     End of DGET39
  *
        END
  