|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 |
SUBROUTINE DSVDCT( N, S, E, SHIFT, NUM )
* * -- LAPACK test routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER N, NUM DOUBLE PRECISION SHIFT * .. * .. Array Arguments .. DOUBLE PRECISION E( * ), S( * ) * .. * * Purpose * ======= * * DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N * tridiagonal matrix T which are less than or equal to SHIFT. T is * formed by putting zeros on the diagonal and making the off-diagonals * equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N). If SHIFT is * positive, NUM is equal to N plus the number of singular values of a * bidiagonal matrix B less than or equal to SHIFT. Here B has diagonal * entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1). * If SHIFT is negative, NUM is equal to the number of singular values * of B greater than or equal to -SHIFT. * * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal * Matrix", Report CS41, Computer Science Dept., Stanford University, * July 21, 1966 * * Arguments * ========= * * N (input) INTEGER * The dimension of the bidiagonal matrix B. * * S (input) DOUBLE PRECISION array, dimension (N) * The diagonal entries of the bidiagonal matrix B. * * E (input) DOUBLE PRECISION array of dimension (N-1) * The superdiagonal entries of the bidiagonal matrix B. * * SHIFT (input) DOUBLE PRECISION * The shift, used as described under Purpose. * * NUM (output) INTEGER * The number of eigenvalues of T less than or equal to SHIFT. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) * .. * .. Local Scalars .. INTEGER I DOUBLE PRECISION M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP, $ TOM, U, UNFL * .. * .. External Functions .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. * .. Executable Statements .. * * Get machine constants * UNFL = 2*DLAMCH( 'Safe minimum' ) OVFL = ONE / UNFL * * Find largest entry * MX = ABS( S( 1 ) ) DO 10 I = 1, N - 1 MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) ) 10 CONTINUE * IF( MX.EQ.ZERO ) THEN IF( SHIFT.LT.ZERO ) THEN NUM = 0 ELSE NUM = 2*N END IF RETURN END IF * * Compute scale factors as in Kahan's report * SUN = SQRT( UNFL ) SSUN = SQRT( SUN ) SOV = SQRT( OVFL ) TOM = SSUN*SOV IF( MX.LE.ONE ) THEN M1 = ONE / MX M2 = TOM ELSE M1 = ONE M2 = TOM / MX END IF * * Begin counting * U = ONE NUM = 0 SSHIFT = ( SHIFT*M1 )*M2 U = -SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF TMP = ( S( 1 )*M1 )*M2 U = -TMP*( TMP / U ) - SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF DO 20 I = 1, N - 1 TMP = ( E( I )*M1 )*M2 U = -TMP*( TMP / U ) - SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF TMP = ( S( I+1 )*M1 )*M2 U = -TMP*( TMP / U ) - SSHIFT IF( U.LE.SUN ) THEN IF( U.LE.ZERO ) THEN NUM = NUM + 1 IF( U.GT.-SUN ) $ U = -SUN ELSE U = SUN END IF END IF 20 CONTINUE RETURN * * End of DSVDCT * END |