![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSTECH 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 DSTECH( N, A, B, EIG, TOL, WORK, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, N 00015 * DOUBLE PRECISION TOL 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION A( * ), B( * ), EIG( * ), WORK( * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> Let T be the tridiagonal matrix with diagonal entries A(1) ,..., 00028 *> A(N) and offdiagonal entries B(1) ,..., B(N-1)). DSTECH checks to 00029 *> see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T. 00030 *> It does this by expanding each EIG(I) into an interval 00031 *> [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if 00032 *> any, and using Sturm sequences to count and verify whether each 00033 *> resulting interval has the correct number of eigenvalues (using 00034 *> DSTECT). Here EPS = TOL*MAZHEPS*MAXEIG, where MACHEPS is the 00035 *> machine precision and MAXEIG is the absolute value of the largest 00036 *> eigenvalue. If each interval contains the correct number of 00037 *> eigenvalues, INFO = 0 is returned, otherwise INFO is the index of 00038 *> the first eigenvalue in the first bad interval. 00039 *> \endverbatim 00040 * 00041 * Arguments: 00042 * ========== 00043 * 00044 *> \param[in] N 00045 *> \verbatim 00046 *> N is INTEGER 00047 *> The dimension of the tridiagonal matrix T. 00048 *> \endverbatim 00049 *> 00050 *> \param[in] A 00051 *> \verbatim 00052 *> A is DOUBLE PRECISION array, dimension (N) 00053 *> The diagonal entries of the tridiagonal matrix T. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] B 00057 *> \verbatim 00058 *> B is DOUBLE PRECISION array, dimension (N-1) 00059 *> The offdiagonal entries of the tridiagonal matrix T. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] EIG 00063 *> \verbatim 00064 *> EIG is DOUBLE PRECISION array, dimension (N) 00065 *> The purported eigenvalues to be checked. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] TOL 00069 *> \verbatim 00070 *> TOL is DOUBLE PRECISION 00071 *> Error tolerance for checking, a multiple of the 00072 *> machine precision. 00073 *> \endverbatim 00074 *> 00075 *> \param[out] WORK 00076 *> \verbatim 00077 *> WORK is DOUBLE PRECISION array, dimension (N) 00078 *> \endverbatim 00079 *> 00080 *> \param[out] INFO 00081 *> \verbatim 00082 *> INFO is INTEGER 00083 *> 0 if the eigenvalues are all correct (to within 00084 *> 1 +- TOL*MAZHEPS*MAXEIG) 00085 *> >0 if the interval containing the INFO-th eigenvalue 00086 *> contains the incorrect number of eigenvalues. 00087 *> \endverbatim 00088 * 00089 * Authors: 00090 * ======== 00091 * 00092 *> \author Univ. of Tennessee 00093 *> \author Univ. of California Berkeley 00094 *> \author Univ. of Colorado Denver 00095 *> \author NAG Ltd. 00096 * 00097 *> \date November 2011 00098 * 00099 *> \ingroup double_eig 00100 * 00101 * ===================================================================== 00102 SUBROUTINE DSTECH( N, A, B, EIG, TOL, WORK, INFO ) 00103 * 00104 * -- LAPACK test routine (version 3.4.0) -- 00105 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00106 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00107 * November 2011 00108 * 00109 * .. Scalar Arguments .. 00110 INTEGER INFO, N 00111 DOUBLE PRECISION TOL 00112 * .. 00113 * .. Array Arguments .. 00114 DOUBLE PRECISION A( * ), B( * ), EIG( * ), WORK( * ) 00115 * .. 00116 * 00117 * ===================================================================== 00118 * 00119 * .. Parameters .. 00120 DOUBLE PRECISION ZERO 00121 PARAMETER ( ZERO = 0.0D0 ) 00122 * .. 00123 * .. Local Scalars .. 00124 INTEGER BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT 00125 DOUBLE PRECISION EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER 00126 * .. 00127 * .. External Functions .. 00128 DOUBLE PRECISION DLAMCH 00129 EXTERNAL DLAMCH 00130 * .. 00131 * .. External Subroutines .. 00132 EXTERNAL DSTECT 00133 * .. 00134 * .. Intrinsic Functions .. 00135 INTRINSIC ABS, MAX 00136 * .. 00137 * .. Executable Statements .. 00138 * 00139 * Check input parameters 00140 * 00141 INFO = 0 00142 IF( N.EQ.0 ) 00143 $ RETURN 00144 IF( N.LT.0 ) THEN 00145 INFO = -1 00146 RETURN 00147 END IF 00148 IF( TOL.LT.ZERO ) THEN 00149 INFO = -5 00150 RETURN 00151 END IF 00152 * 00153 * Get machine constants 00154 * 00155 EPS = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00156 UNFLEP = DLAMCH( 'Safe minimum' ) / EPS 00157 EPS = TOL*EPS 00158 * 00159 * Compute maximum absolute eigenvalue, error tolerance 00160 * 00161 MX = ABS( EIG( 1 ) ) 00162 DO 10 I = 2, N 00163 MX = MAX( MX, ABS( EIG( I ) ) ) 00164 10 CONTINUE 00165 EPS = MAX( EPS*MX, UNFLEP ) 00166 * 00167 * Sort eigenvalues from EIG into WORK 00168 * 00169 DO 20 I = 1, N 00170 WORK( I ) = EIG( I ) 00171 20 CONTINUE 00172 DO 40 I = 1, N - 1 00173 ISUB = 1 00174 EMIN = WORK( 1 ) 00175 DO 30 J = 2, N + 1 - I 00176 IF( WORK( J ).LT.EMIN ) THEN 00177 ISUB = J 00178 EMIN = WORK( J ) 00179 END IF 00180 30 CONTINUE 00181 IF( ISUB.NE.N+1-I ) THEN 00182 WORK( ISUB ) = WORK( N+1-I ) 00183 WORK( N+1-I ) = EMIN 00184 END IF 00185 40 CONTINUE 00186 * 00187 * TPNT points to singular value at right endpoint of interval 00188 * BPNT points to singular value at left endpoint of interval 00189 * 00190 TPNT = 1 00191 BPNT = 1 00192 * 00193 * Begin loop over all intervals 00194 * 00195 50 CONTINUE 00196 UPPER = WORK( TPNT ) + EPS 00197 LOWER = WORK( BPNT ) - EPS 00198 * 00199 * Begin loop merging overlapping intervals 00200 * 00201 60 CONTINUE 00202 IF( BPNT.EQ.N ) 00203 $ GO TO 70 00204 TUPPR = WORK( BPNT+1 ) + EPS 00205 IF( TUPPR.LT.LOWER ) 00206 $ GO TO 70 00207 * 00208 * Merge 00209 * 00210 BPNT = BPNT + 1 00211 LOWER = WORK( BPNT ) - EPS 00212 GO TO 60 00213 70 CONTINUE 00214 * 00215 * Count singular values in interval [ LOWER, UPPER ] 00216 * 00217 CALL DSTECT( N, A, B, LOWER, NUML ) 00218 CALL DSTECT( N, A, B, UPPER, NUMU ) 00219 COUNT = NUMU - NUML 00220 IF( COUNT.NE.BPNT-TPNT+1 ) THEN 00221 * 00222 * Wrong number of singular values in interval 00223 * 00224 INFO = TPNT 00225 GO TO 80 00226 END IF 00227 TPNT = BPNT + 1 00228 BPNT = TPNT 00229 IF( TPNT.LE.N ) 00230 $ GO TO 50 00231 80 CONTINUE 00232 RETURN 00233 * 00234 * End of DSTECH 00235 * 00236 END