00001 SUBROUTINE DLARRK( N, IW, GL, GU, 00002 $ D, E2, PIVMIN, RELTOL, W, WERR, INFO) 00003 IMPLICIT NONE 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER INFO, IW, N 00012 DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION D( * ), E2( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLARRK computes one eigenvalue of a symmetric tridiagonal 00022 * matrix T to suitable accuracy. This is an auxiliary code to be 00023 * called from DSTEMR. 00024 * 00025 * To avoid overflow, the matrix must be scaled so that its 00026 * largest element is no greater than overflow**(1/2) * 00027 * underflow**(1/4) in absolute value, and for greatest 00028 * accuracy, it should not be much smaller than that. 00029 * 00030 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal 00031 * Matrix", Report CS41, Computer Science Dept., Stanford 00032 * University, July 21, 1966. 00033 * 00034 * Arguments 00035 * ========= 00036 * 00037 * N (input) INTEGER 00038 * The order of the tridiagonal matrix T. N >= 0. 00039 * 00040 * IW (input) INTEGER 00041 * The index of the eigenvalues to be returned. 00042 * 00043 * GL (input) DOUBLE PRECISION 00044 * GU (input) DOUBLE PRECISION 00045 * An upper and a lower bound on the eigenvalue. 00046 * 00047 * D (input) DOUBLE PRECISION array, dimension (N) 00048 * The n diagonal elements of the tridiagonal matrix T. 00049 * 00050 * E2 (input) DOUBLE PRECISION array, dimension (N-1) 00051 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T. 00052 * 00053 * PIVMIN (input) DOUBLE PRECISION 00054 * The minimum pivot allowed in the Sturm sequence for T. 00055 * 00056 * RELTOL (input) DOUBLE PRECISION 00057 * The minimum relative width of an interval. When an interval 00058 * is narrower than RELTOL times the larger (in 00059 * magnitude) endpoint, then it is considered to be 00060 * sufficiently small, i.e., converged. Note: this should 00061 * always be at least radix*machine epsilon. 00062 * 00063 * W (output) DOUBLE PRECISION 00064 * 00065 * WERR (output) DOUBLE PRECISION 00066 * The error bound on the corresponding eigenvalue approximation 00067 * in W. 00068 * 00069 * INFO (output) INTEGER 00070 * = 0: Eigenvalue converged 00071 * = -1: Eigenvalue did NOT converge 00072 * 00073 * Internal Parameters 00074 * =================== 00075 * 00076 * FUDGE DOUBLE PRECISION, default = 2 00077 * A "fudge factor" to widen the Gershgorin intervals. 00078 * 00079 * ===================================================================== 00080 * 00081 * .. Parameters .. 00082 DOUBLE PRECISION FUDGE, HALF, TWO, ZERO 00083 PARAMETER ( HALF = 0.5D0, TWO = 2.0D0, 00084 $ FUDGE = TWO, ZERO = 0.0D0 ) 00085 * .. 00086 * .. Local Scalars .. 00087 INTEGER I, IT, ITMAX, NEGCNT 00088 DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1, 00089 $ TMP2, TNORM 00090 * .. 00091 * .. External Functions .. 00092 DOUBLE PRECISION DLAMCH 00093 EXTERNAL DLAMCH 00094 * .. 00095 * .. Intrinsic Functions .. 00096 INTRINSIC ABS, INT, LOG, MAX 00097 * .. 00098 * .. Executable Statements .. 00099 * 00100 * Get machine constants 00101 EPS = DLAMCH( 'P' ) 00102 00103 TNORM = MAX( ABS( GL ), ABS( GU ) ) 00104 RTOLI = RELTOL 00105 ATOLI = FUDGE*TWO*PIVMIN 00106 00107 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / 00108 $ LOG( TWO ) ) + 2 00109 00110 INFO = -1 00111 00112 LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN 00113 RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN 00114 IT = 0 00115 00116 10 CONTINUE 00117 * 00118 * Check if interval converged or maximum number of iterations reached 00119 * 00120 TMP1 = ABS( RIGHT - LEFT ) 00121 TMP2 = MAX( ABS(RIGHT), ABS(LEFT) ) 00122 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN 00123 INFO = 0 00124 GOTO 30 00125 ENDIF 00126 IF(IT.GT.ITMAX) 00127 $ GOTO 30 00128 00129 * 00130 * Count number of negative pivots for mid-point 00131 * 00132 IT = IT + 1 00133 MID = HALF * (LEFT + RIGHT) 00134 NEGCNT = 0 00135 TMP1 = D( 1 ) - MID 00136 IF( ABS( TMP1 ).LT.PIVMIN ) 00137 $ TMP1 = -PIVMIN 00138 IF( TMP1.LE.ZERO ) 00139 $ NEGCNT = NEGCNT + 1 00140 * 00141 DO 20 I = 2, N 00142 TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID 00143 IF( ABS( TMP1 ).LT.PIVMIN ) 00144 $ TMP1 = -PIVMIN 00145 IF( TMP1.LE.ZERO ) 00146 $ NEGCNT = NEGCNT + 1 00147 20 CONTINUE 00148 00149 IF(NEGCNT.GE.IW) THEN 00150 RIGHT = MID 00151 ELSE 00152 LEFT = MID 00153 ENDIF 00154 GOTO 10 00155 00156 30 CONTINUE 00157 * 00158 * Converged or maximum number of iterations reached 00159 * 00160 W = HALF * (LEFT + RIGHT) 00161 WERR = HALF * ABS( RIGHT - LEFT ) 00162 00163 RETURN 00164 * 00165 * End of DLARRK 00166 * 00167 END
1.7.2