![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSVDCT 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DSVDCT( N, S, E, SHIFT, NUM ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER N, NUM 00015 * DOUBLE PRECISION SHIFT 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION E( * ), S( * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DSVDCT counts the number NUM of eigenvalues of a 2*N by 2*N 00028 *> tridiagonal matrix T which are less than or equal to SHIFT. T is 00029 *> formed by putting zeros on the diagonal and making the off-diagonals 00030 *> equal to S(1), E(1), S(2), E(2), ... , E(N-1), S(N). If SHIFT is 00031 *> positive, NUM is equal to N plus the number of singular values of a 00032 *> bidiagonal matrix B less than or equal to SHIFT. Here B has diagonal 00033 *> entries S(1), ..., S(N) and superdiagonal entries E(1), ... E(N-1). 00034 *> If SHIFT is negative, NUM is equal to the number of singular values 00035 *> of B greater than or equal to -SHIFT. 00036 *> 00037 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00038 *> Matrix", Report CS41, Computer Science Dept., Stanford University, 00039 *> July 21, 1966 00040 *> \endverbatim 00041 * 00042 * Arguments: 00043 * ========== 00044 * 00045 *> \param[in] N 00046 *> \verbatim 00047 *> N is INTEGER 00048 *> The dimension of the bidiagonal matrix B. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] S 00052 *> \verbatim 00053 *> S is DOUBLE PRECISION array, dimension (N) 00054 *> The diagonal entries of the bidiagonal matrix B. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] E 00058 *> \verbatim 00059 *> E is DOUBLE PRECISION array of dimension (N-1) 00060 *> The superdiagonal entries of the bidiagonal matrix B. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] SHIFT 00064 *> \verbatim 00065 *> SHIFT is DOUBLE PRECISION 00066 *> The shift, used as described under Purpose. 00067 *> \endverbatim 00068 *> 00069 *> \param[out] NUM 00070 *> \verbatim 00071 *> NUM is INTEGER 00072 *> The number of eigenvalues of T less than or equal to SHIFT. 00073 *> \endverbatim 00074 * 00075 * Authors: 00076 * ======== 00077 * 00078 *> \author Univ. of Tennessee 00079 *> \author Univ. of California Berkeley 00080 *> \author Univ. of Colorado Denver 00081 *> \author NAG Ltd. 00082 * 00083 *> \date November 2011 00084 * 00085 *> \ingroup double_eig 00086 * 00087 * ===================================================================== 00088 SUBROUTINE DSVDCT( N, S, E, SHIFT, NUM ) 00089 * 00090 * -- LAPACK test routine (version 3.4.0) -- 00091 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00092 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00093 * November 2011 00094 * 00095 * .. Scalar Arguments .. 00096 INTEGER N, NUM 00097 DOUBLE PRECISION SHIFT 00098 * .. 00099 * .. Array Arguments .. 00100 DOUBLE PRECISION E( * ), S( * ) 00101 * .. 00102 * 00103 * ===================================================================== 00104 * 00105 * .. Parameters .. 00106 DOUBLE PRECISION ONE 00107 PARAMETER ( ONE = 1.0D0 ) 00108 DOUBLE PRECISION ZERO 00109 PARAMETER ( ZERO = 0.0D0 ) 00110 * .. 00111 * .. Local Scalars .. 00112 INTEGER I 00113 DOUBLE PRECISION M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP, 00114 $ TOM, U, UNFL 00115 * .. 00116 * .. External Functions .. 00117 DOUBLE PRECISION DLAMCH 00118 EXTERNAL DLAMCH 00119 * .. 00120 * .. Intrinsic Functions .. 00121 INTRINSIC ABS, MAX, SQRT 00122 * .. 00123 * .. Executable Statements .. 00124 * 00125 * Get machine constants 00126 * 00127 UNFL = 2*DLAMCH( 'Safe minimum' ) 00128 OVFL = ONE / UNFL 00129 * 00130 * Find largest entry 00131 * 00132 MX = ABS( S( 1 ) ) 00133 DO 10 I = 1, N - 1 00134 MX = MAX( MX, ABS( S( I+1 ) ), ABS( E( I ) ) ) 00135 10 CONTINUE 00136 * 00137 IF( MX.EQ.ZERO ) THEN 00138 IF( SHIFT.LT.ZERO ) THEN 00139 NUM = 0 00140 ELSE 00141 NUM = 2*N 00142 END IF 00143 RETURN 00144 END IF 00145 * 00146 * Compute scale factors as in Kahan's report 00147 * 00148 SUN = SQRT( UNFL ) 00149 SSUN = SQRT( SUN ) 00150 SOV = SQRT( OVFL ) 00151 TOM = SSUN*SOV 00152 IF( MX.LE.ONE ) THEN 00153 M1 = ONE / MX 00154 M2 = TOM 00155 ELSE 00156 M1 = ONE 00157 M2 = TOM / MX 00158 END IF 00159 * 00160 * Begin counting 00161 * 00162 U = ONE 00163 NUM = 0 00164 SSHIFT = ( SHIFT*M1 )*M2 00165 U = -SSHIFT 00166 IF( U.LE.SUN ) THEN 00167 IF( U.LE.ZERO ) THEN 00168 NUM = NUM + 1 00169 IF( U.GT.-SUN ) 00170 $ U = -SUN 00171 ELSE 00172 U = SUN 00173 END IF 00174 END IF 00175 TMP = ( S( 1 )*M1 )*M2 00176 U = -TMP*( TMP / U ) - SSHIFT 00177 IF( U.LE.SUN ) THEN 00178 IF( U.LE.ZERO ) THEN 00179 NUM = NUM + 1 00180 IF( U.GT.-SUN ) 00181 $ U = -SUN 00182 ELSE 00183 U = SUN 00184 END IF 00185 END IF 00186 DO 20 I = 1, N - 1 00187 TMP = ( E( I )*M1 )*M2 00188 U = -TMP*( TMP / U ) - SSHIFT 00189 IF( U.LE.SUN ) THEN 00190 IF( U.LE.ZERO ) THEN 00191 NUM = NUM + 1 00192 IF( U.GT.-SUN ) 00193 $ U = -SUN 00194 ELSE 00195 U = SUN 00196 END IF 00197 END IF 00198 TMP = ( S( I+1 )*M1 )*M2 00199 U = -TMP*( TMP / U ) - SSHIFT 00200 IF( U.LE.SUN ) THEN 00201 IF( U.LE.ZERO ) THEN 00202 NUM = NUM + 1 00203 IF( U.GT.-SUN ) 00204 $ U = -SUN 00205 ELSE 00206 U = SUN 00207 END IF 00208 END IF 00209 20 CONTINUE 00210 RETURN 00211 * 00212 * End of DSVDCT 00213 * 00214 END