TESTING/EIG/dstech.f

   1 *> \brief \b DSTECH
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            INFO, N
  15 *       DOUBLE PRECISION   TOL
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       DOUBLE PRECISION   A( * ), B( * ), EIG( * ), WORK( * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *>    Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
  28 *>    A(N) and offdiagonal entries B(1) ,..., B(N-1)).  DSTECH checks to
  29 *>    see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
  30 *>    It does this by expanding each EIG(I) into an interval
  31 *>    [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
  32 *>    any, and using Sturm sequences to count and verify whether each
  33 *>    resulting interval has the correct number of eigenvalues (using
  34 *>    DSTECT).  Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the
  35 *>    machine precision and MAXEIG is the absolute value of the largest
  36 *>    eigenvalue. If each interval contains the correct number of
  37 *>    eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
  38 *>    the first eigenvalue in the first bad interval.
  39 *> \endverbatim
  40 *
  41 *  Arguments:
  42 *  ==========
  43 *
  44 *> \param[in] N
  45 *> \verbatim
  46 *>          N is INTEGER
  47 *>          The dimension of the tridiagonal matrix T.
  48 *> \endverbatim
  49 *>
  50 *> \param[in] A
  51 *> \verbatim
  52 *>          A is DOUBLE PRECISION array, dimension (N)
  53 *>          The diagonal entries of the tridiagonal matrix T.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] B
  57 *> \verbatim
  58 *>          B is DOUBLE PRECISION array, dimension (N-1)
  59 *>          The offdiagonal entries of the tridiagonal matrix T.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] EIG
  63 *> \verbatim
  64 *>          EIG is DOUBLE PRECISION array, dimension (N)
  65 *>          The purported eigenvalues to be checked.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] TOL
  69 *> \verbatim
  70 *>          TOL is DOUBLE PRECISION
  71 *>          Error tolerance for checking, a multiple of the
  72 *>          machine precision.
  73 *> \endverbatim
  74 *>
  75 *> \param[out] WORK
  76 *> \verbatim
  77 *>          WORK is DOUBLE PRECISION array, dimension (N)
  78 *> \endverbatim
  79 *>
  80 *> \param[out] INFO
  81 *> \verbatim
  82 *>          INFO is INTEGER
  83 *>          0  if the eigenvalues are all correct (to within
  84 *>             1 +- TOL*MAZHEPS*MAXEIG)
  85 *>          >0 if the interval containing the INFO-th eigenvalue
  86 *>             contains the incorrect number of eigenvalues.
  87 *> \endverbatim
  88 *
  89 *  Authors:
  90 *  ========
  91 *
  92 *> \author Univ. of Tennessee
  93 *> \author Univ. of California Berkeley
  94 *> \author Univ. of Colorado Denver
  95 *> \author NAG Ltd.
  96 *
  97 *> \date November 2011
  98 *
  99 *> \ingroup double_eig
 100 *
 101 *  =====================================================================
 102       SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO )
 103 *
 104 *  -- LAPACK test routine (version 3.4.0) --
 105 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 106 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 107 *     November 2011
 108 *
 109 *     .. Scalar Arguments ..
 110       INTEGER            INFO, N
 111       DOUBLE PRECISION   TOL
 112 *     ..
 113 *     .. Array Arguments ..
 114       DOUBLE PRECISION   A( * ), B( * ), EIG( * ), WORK( * )
 115 *     ..
 116 *
 117 *  =====================================================================
 118 *
 119 *     .. Parameters ..
 120       DOUBLE PRECISION   ZERO
 121       PARAMETER          ( ZERO = 0.0D0 )
 122 *     ..
 123 *     .. Local Scalars ..
 124       INTEGER            BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT
 125       DOUBLE PRECISION   EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER
 126 *     ..
 127 *     .. External Functions ..
 128       DOUBLE PRECISION   DLAMCH
 129       EXTERNAL           DLAMCH
 130 *     ..
 131 *     .. External Subroutines ..
 132       EXTERNAL           DSTECT
 133 *     ..
 134 *     .. Intrinsic Functions ..
 135       INTRINSIC          ABS, MAX
 136 *     ..
 137 *     .. Executable Statements ..
 138 *
 139 *     Check input parameters
 140 *
 141       INFO = 0
 142       IF( N.EQ.0 )
 143      $   RETURN
 144       IF( N.LT.0 ) THEN
 145          INFO = -1
 146          RETURN
 147       END IF
 148       IF( TOL.LT.ZERO ) THEN
 149          INFO = -5
 150          RETURN
 151       END IF
 152 *
 153 *     Get machine constants
 154 *
 155       EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
 156       UNFLEP = DLAMCH( 'Safe minimum' ) / EPS
 157       EPS = TOL*EPS
 158 *
 159 *     Compute maximum absolute eigenvalue, error tolerance
 160 *
 161       MX = ABS( EIG( 1 ) )
 162       DO 10 I = 2, N
 163          MX = MAX( MX, ABS( EIG( I ) ) )
 164    10 CONTINUE
 165       EPS = MAX( EPS*MX, UNFLEP )
 166 *
 167 *     Sort eigenvalues from EIG into WORK
 168 *
 169       DO 20 I = 1, N
 170          WORK( I ) = EIG( I )
 171    20 CONTINUE
 172       DO 40 I = 1, N - 1
 173          ISUB = 1
 174          EMIN = WORK( 1 )
 175          DO 30 J = 2, N + 1 - I
 176             IF( WORK( J ).LT.EMIN ) THEN
 177                ISUB = J
 178                EMIN = WORK( J )
 179             END IF
 180    30    CONTINUE
 181          IF( ISUB.NE.N+1-I ) THEN
 182             WORK( ISUB ) = WORK( N+1-I )
 183             WORK( N+1-I ) = EMIN
 184          END IF
 185    40 CONTINUE
 186 *
 187 *     TPNT points to singular value at right endpoint of interval
 188 *     BPNT points to singular value at left  endpoint of interval
 189 *
 190       TPNT = 1
 191       BPNT = 1
 192 *
 193 *     Begin loop over all intervals
 194 *
 195    50 CONTINUE
 196       UPPER = WORK( TPNT ) + EPS
 197       LOWER = WORK( BPNT ) - EPS
 198 *
 199 *     Begin loop merging overlapping intervals
 200 *
 201    60 CONTINUE
 202       IF( BPNT.EQ.N )
 203      $   GO TO 70
 204       TUPPR = WORK( BPNT+1 ) + EPS
 205       IF( TUPPR.LT.LOWER )
 206      $   GO TO 70
 207 *
 208 *     Merge
 209 *
 210       BPNT = BPNT + 1
 211       LOWER = WORK( BPNT ) - EPS
 212       GO TO 60
 213    70 CONTINUE
 214 *
 215 *     Count singular values in interval [ LOWER, UPPER ]
 216 *
 217       CALL DSTECT( N, A, B, LOWER, NUML )
 218       CALL DSTECT( N, A, B, UPPER, NUMU )
 219       COUNT = NUMU - NUML
 220       IF( COUNT.NE.BPNT-TPNT+1 ) THEN
 221 *
 222 *        Wrong number of singular values in interval
 223 *
 224          INFO = TPNT
 225          GO TO 80
 226       END IF
 227       TPNT = BPNT + 1
 228       BPNT = TPNT
 229       IF( TPNT.LE.N )
 230      $   GO TO 50
 231    80 CONTINUE
 232       RETURN
 233 *
 234 *     End of DSTECH
 235 *
 236       END