![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTECT 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 DSTECT( N, A, B, SHIFT, NUM ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER N, NUM 00015 * DOUBLE PRECISION SHIFT 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION A( * ), B( * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DSTECT counts the number NUM of eigenvalues of a tridiagonal 00028 *> matrix T which are less than or equal to SHIFT. T has 00029 *> diagonal entries A(1), ... , A(N), and offdiagonal entries 00030 *> B(1), ..., B(N-1). 00031 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00032 *> Matrix", Report CS41, Computer Science Dept., Stanford 00033 *> University, July 21, 1966 00034 *> \endverbatim 00035 * 00036 * Arguments: 00037 * ========== 00038 * 00039 *> \param[in] N 00040 *> \verbatim 00041 *> N is INTEGER 00042 *> The dimension of the tridiagonal matrix T. 00043 *> \endverbatim 00044 *> 00045 *> \param[in] A 00046 *> \verbatim 00047 *> A is DOUBLE PRECISION array, dimension (N) 00048 *> The diagonal entries of the tridiagonal matrix T. 00049 *> \endverbatim 00050 *> 00051 *> \param[in] B 00052 *> \verbatim 00053 *> B is DOUBLE PRECISION array, dimension (N-1) 00054 *> The offdiagonal entries of the tridiagonal matrix T. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] SHIFT 00058 *> \verbatim 00059 *> SHIFT is DOUBLE PRECISION 00060 *> The shift, used as described under Purpose. 00061 *> \endverbatim 00062 *> 00063 *> \param[out] NUM 00064 *> \verbatim 00065 *> NUM is INTEGER 00066 *> The number of eigenvalues of T less than or equal 00067 *> to SHIFT. 00068 *> \endverbatim 00069 * 00070 * Authors: 00071 * ======== 00072 * 00073 *> \author Univ. of Tennessee 00074 *> \author Univ. of California Berkeley 00075 *> \author Univ. of Colorado Denver 00076 *> \author NAG Ltd. 00077 * 00078 *> \date November 2011 00079 * 00080 *> \ingroup double_eig 00081 * 00082 * ===================================================================== 00083 SUBROUTINE DSTECT( N, A, B, SHIFT, NUM ) 00084 * 00085 * -- LAPACK test routine (version 3.4.0) -- 00086 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00087 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00088 * November 2011 00089 * 00090 * .. Scalar Arguments .. 00091 INTEGER N, NUM 00092 DOUBLE PRECISION SHIFT 00093 * .. 00094 * .. Array Arguments .. 00095 DOUBLE PRECISION A( * ), B( * ) 00096 * .. 00097 * 00098 * ===================================================================== 00099 * 00100 * .. Parameters .. 00101 DOUBLE PRECISION ZERO, ONE, THREE 00102 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, THREE = 3.0D0 ) 00103 * .. 00104 * .. Local Scalars .. 00105 INTEGER I 00106 DOUBLE PRECISION M1, M2, MX, OVFL, SOV, SSHIFT, SSUN, SUN, TMP, 00107 $ TOM, U, UNFL 00108 * .. 00109 * .. External Functions .. 00110 DOUBLE PRECISION DLAMCH 00111 EXTERNAL DLAMCH 00112 * .. 00113 * .. Intrinsic Functions .. 00114 INTRINSIC ABS, MAX, SQRT 00115 * .. 00116 * .. Executable Statements .. 00117 * 00118 * Get machine constants 00119 * 00120 UNFL = DLAMCH( 'Safe minimum' ) 00121 OVFL = DLAMCH( 'Overflow' ) 00122 * 00123 * Find largest entry 00124 * 00125 MX = ABS( A( 1 ) ) 00126 DO 10 I = 1, N - 1 00127 MX = MAX( MX, ABS( A( I+1 ) ), ABS( B( I ) ) ) 00128 10 CONTINUE 00129 * 00130 * Handle easy cases, including zero matrix 00131 * 00132 IF( SHIFT.GE.THREE*MX ) THEN 00133 NUM = N 00134 RETURN 00135 END IF 00136 IF( SHIFT.LT.-THREE*MX ) THEN 00137 NUM = 0 00138 RETURN 00139 END IF 00140 * 00141 * Compute scale factors as in Kahan's report 00142 * At this point, MX .NE. 0 so we can divide by it 00143 * 00144 SUN = SQRT( UNFL ) 00145 SSUN = SQRT( SUN ) 00146 SOV = SQRT( OVFL ) 00147 TOM = SSUN*SOV 00148 IF( MX.LE.ONE ) THEN 00149 M1 = ONE / MX 00150 M2 = TOM 00151 ELSE 00152 M1 = ONE 00153 M2 = TOM / MX 00154 END IF 00155 * 00156 * Begin counting 00157 * 00158 NUM = 0 00159 SSHIFT = ( SHIFT*M1 )*M2 00160 U = ( A( 1 )*M1 )*M2 - SSHIFT 00161 IF( U.LE.SUN ) THEN 00162 IF( U.LE.ZERO ) THEN 00163 NUM = NUM + 1 00164 IF( U.GT.-SUN ) 00165 $ U = -SUN 00166 ELSE 00167 U = SUN 00168 END IF 00169 END IF 00170 DO 20 I = 2, N 00171 TMP = ( B( I-1 )*M1 )*M2 00172 U = ( ( A( I )*M1 )*M2-TMP*( TMP / U ) ) - SSHIFT 00173 IF( U.LE.SUN ) THEN 00174 IF( U.LE.ZERO ) THEN 00175 NUM = NUM + 1 00176 IF( U.GT.-SUN ) 00177 $ U = -SUN 00178 ELSE 00179 U = SUN 00180 END IF 00181 END IF 00182 20 CONTINUE 00183 RETURN 00184 * 00185 * End of DSTECT 00186 * 00187 END