TESTING/EIG/dsvdct.f

   1 *> \brief \b DSVDCT
   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 DSVDCT( N, S, E, SHIFT, NUM )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            N, NUM
  15 *       DOUBLE PRECISION   SHIFT
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       DOUBLE PRECISION   E( * ), S( * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *> DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N
  28 *> tridiagonal matrix T which are less than or equal to SHIFT.  T is
  29 *> formed by putting zeros on the diagonal and making the off-diagonals
  30 *> equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N).  If SHIFT is
  31 *> positive, NUM is equal to N plus the number of singular values of a
  32 *> bidiagonal matrix B less than or equal to SHIFT.  Here B has diagonal
  33 *> entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1).
  34 *> If SHIFT is negative, NUM is equal to the number of singular values
  35 *> of B greater than or equal to -SHIFT.
  36 *>
  37 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
  38 *> Matrix", Report CS41, Computer Science Dept., Stanford University,
  39 *> July 21, 1966
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] N
  46 *> \verbatim
  47 *>          N is INTEGER
  48 *>          The dimension of the bidiagonal matrix B.
  49 *> \endverbatim
  50 *>
  51 *> \param[in] S
  52 *> \verbatim
  53 *>          S is DOUBLE PRECISION array, dimension (N)
  54 *>          The diagonal entries of the bidiagonal matrix B.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] E
  58 *> \verbatim
  59 *>          E is DOUBLE PRECISION array of dimension (N-1)
  60 *>          The superdiagonal entries of the bidiagonal matrix B.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] SHIFT
  64 *> \verbatim
  65 *>          SHIFT is DOUBLE PRECISION
  66 *>          The shift, used as described under Purpose.
  67 *> \endverbatim
  68 *>
  69 *> \param[out] NUM
  70 *> \verbatim
  71 *>          NUM is INTEGER
  72 *>          The number of eigenvalues of T less than or equal to SHIFT.
  73 *> \endverbatim
  74 *
  75 *  Authors:
  76 *  ========
  77 *
  78 *> \author Univ. of Tennessee
  79 *> \author Univ. of California Berkeley
  80 *> \author Univ. of Colorado Denver
  81 *> \author NAG Ltd.
  82 *
  83 *> \date November 2011
  84 *
  85 *> \ingroup double_eig
  86 *
  87 *  =====================================================================
  88       SUBROUTINE DSVDCT( N, S, E, SHIFT, NUM )
  89 *
  90 *  -- LAPACK test routine (version 3.4.0) --
  91 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  92 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  93 *     November 2011
  94 *
  95 *     .. Scalar Arguments ..
  96       INTEGER            N, NUM
  97       DOUBLE PRECISION   SHIFT
  98 *     ..
  99 *     .. Array Arguments ..
 100       DOUBLE PRECISION   E( * ), S( * )
 101 *     ..
 102 *
 103 *  =====================================================================
 104 *
 105 *     .. Parameters ..
 106       DOUBLE PRECISION   ONE
 107       PARAMETER          ( ONE = 1.0D0 )
 108       DOUBLE PRECISION   ZERO
 109       PARAMETER          ( ZERO = 0.0D0 )
 110 *     ..
 111 *     .. Local Scalars ..
 112       INTEGER            I
 113       DOUBLE PRECISION   M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP,
 114      $                   TOM, U, UNFL
 115 *     ..
 116 *     .. External Functions ..
 117       DOUBLE PRECISION   DLAMCH
 118       EXTERNAL           DLAMCH
 119 *     ..
 120 *     .. Intrinsic Functions ..
 121       INTRINSIC          ABS, MAX, SQRT
 122 *     ..
 123 *     .. Executable Statements ..
 124 *
 125 *     Get machine constants
 126 *
 127       UNFL = 2*DLAMCH( 'Safe minimum' )
 128       OVFL = ONE / UNFL
 129 *
 130 *     Find largest entry
 131 *
 132       MX = ABS( S( 1 ) )
 133       DO 10 I = 1, N - 1
 134          MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) )
 135    10 CONTINUE
 136 *
 137       IF( MX.EQ.ZERO ) THEN
 138          IF( SHIFT.LT.ZERO ) THEN
 139             NUM = 0
 140          ELSE
 141             NUM = 2*N
 142          END IF
 143          RETURN
 144       END IF
 145 *
 146 *     Compute scale factors as in Kahan's report
 147 *
 148       SUN = SQRT( UNFL )
 149       SSUN = SQRT( SUN )
 150       SOV = SQRT( OVFL )
 151       TOM = SSUN*SOV
 152       IF( MX.LE.ONE ) THEN
 153          M1 = ONE / MX
 154          M2 = TOM
 155       ELSE
 156          M1 = ONE
 157          M2 = TOM / MX
 158       END IF
 159 *
 160 *     Begin counting
 161 *
 162       U = ONE
 163       NUM = 0
 164       SSHIFT = ( SHIFT*M1 )*M2
 165       U = -SSHIFT
 166       IF( U.LE.SUN ) THEN
 167          IF( U.LE.ZERO ) THEN
 168             NUM = NUM + 1
 169             IF( U.GT.-SUN )
 170      $         U = -SUN
 171          ELSE
 172             U = SUN
 173          END IF
 174       END IF
 175       TMP = ( S( 1 )*M1 )*M2
 176       U = -TMP*( TMP / U ) - SSHIFT
 177       IF( U.LE.SUN ) THEN
 178          IF( U.LE.ZERO ) THEN
 179             NUM = NUM + 1
 180             IF( U.GT.-SUN )
 181      $         U = -SUN
 182          ELSE
 183             U = SUN
 184          END IF
 185       END IF
 186       DO 20 I = 1, N - 1
 187          TMP = ( E( I )*M1 )*M2
 188          U = -TMP*( TMP / U ) - SSHIFT
 189          IF( U.LE.SUN ) THEN
 190             IF( U.LE.ZERO ) THEN
 191                NUM = NUM + 1
 192                IF( U.GT.-SUN )
 193      $            U = -SUN
 194             ELSE
 195                U = SUN
 196             END IF
 197          END IF
 198          TMP = ( S( I+1 )*M1 )*M2
 199          U = -TMP*( TMP / U ) - SSHIFT
 200          IF( U.LE.SUN ) THEN
 201             IF( U.LE.ZERO ) THEN
 202                NUM = NUM + 1
 203                IF( U.GT.-SUN )
 204      $            U = -SUN
 205             ELSE
 206                U = SUN
 207             END IF
 208          END IF
 209    20 CONTINUE
 210       RETURN
 211 *
 212 *     End of DSVDCT
 213 *
 214       END