OpenBLAS/lapack-netlib/TESTING/EIG/dget39.f

*> \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 )
*
*     Initalization
*
      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