1 SUBROUTINE DLARRB( N, D, L, LD, LLD, IFIRST, ILAST, SIGMA, RELTOL, 2 $ W, WGAP, WERR, WORK, IWORK, INFO ) 3* 4* -- LAPACK auxiliary routine (instru to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER IFIRST, ILAST, INFO, N 11 DOUBLE PRECISION RELTOL, SIGMA 12* .. 13* .. Array Arguments .. 14 INTEGER IWORK( * ) 15 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ), W( * ), 16 $ WERR( * ), WGAP( * ), WORK( * ) 17* .. 18* Common block to return operation count 19* .. Common blocks .. 20 COMMON / LATIME / OPS, ITCNT 21* .. 22* .. Scalars in Common .. 23 DOUBLE PRECISION ITCNT, OPS 24* .. 25* 26* Purpose 27* ======= 28* 29* Given the relatively robust representation(RRR) L D L^T, DLARRB 30* does ``limited'' bisection to locate the eigenvalues of L D L^T, 31* W( IFIRST ) thru' W( ILAST ), to more accuracy. Intervals 32* [left, right] are maintained by storing their mid-points and 33* semi-widths in the arrays W and WERR respectively. 34* 35* Arguments 36* ========= 37* 38* N (input) INTEGER 39* The order of the matrix. 40* 41* D (input) DOUBLE PRECISION array, dimension (N) 42* The n diagonal elements of the diagonal matrix D. 43* 44* L (input) DOUBLE PRECISION array, dimension (N-1) 45* The n-1 subdiagonal elements of the unit bidiagonal matrix L. 46* 47* LD (input) DOUBLE PRECISION array, dimension (N-1) 48* The n-1 elements L(i)*D(i). 49* 50* LLD (input) DOUBLE PRECISION array, dimension (N-1) 51* The n-1 elements L(i)*L(i)*D(i). 52* 53* IFIRST (input) INTEGER 54* The index of the first eigenvalue in the cluster. 55* 56* ILAST (input) INTEGER 57* The index of the last eigenvalue in the cluster. 58* 59* SIGMA (input) DOUBLE PRECISION 60* The shift used to form L D L^T (see DLARRF). 61* 62* RELTOL (input) DOUBLE PRECISION 63* The relative tolerance. 64* 65* W (input/output) DOUBLE PRECISION array, dimension (N) 66* On input, W( IFIRST ) thru' W( ILAST ) are estimates of the 67* corresponding eigenvalues of L D L^T. 68* On output, these estimates are ``refined''. 69* 70* WGAP (input/output) DOUBLE PRECISION array, dimension (N) 71* The gaps between the eigenvalues of L D L^T. Very small 72* gaps are changed on output. 73* 74* WERR (input/output) DOUBLE PRECISION array, dimension (N) 75* On input, WERR( IFIRST ) thru' WERR( ILAST ) are the errors 76* in the estimates W( IFIRST ) thru' W( ILAST ). 77* On output, these are the ``refined'' errors. 78* 79*****Reminder to Inder --- WORK is never used in this subroutine ***** 80* WORK (input) DOUBLE PRECISION array, dimension (???) 81* Workspace. 82* 83* IWORK (input) INTEGER array, dimension (2*N) 84* Workspace. 85* 86*****Reminder to Inder --- INFO is never set in this subroutine ****** 87* INFO (output) INTEGER 88* Error flag. 89* 90* Further Details 91* =============== 92* 93* Based on contributions by 94* Inderjit Dhillon, IBM Almaden, USA 95* Osni Marques, LBNL/NERSC, USA 96* 97* ===================================================================== 98* 99* .. Parameters .. 100 DOUBLE PRECISION ZERO, TWO, HALF 101 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, HALF = 0.5D0 ) 102* .. 103* .. Local Scalars .. 104 INTEGER CNT, I, I1, I2, INITI1, INITI2, J, K, NCNVRG, 105 $ NEIG, NINT, NRIGHT, OLNINT 106 DOUBLE PRECISION DELTA, EPS, GAP, LEFT, MID, PERT, RIGHT, S, 107 $ THRESH, TMP, WIDTH 108* .. 109* .. External Functions .. 110 DOUBLE PRECISION DLAMCH 111 EXTERNAL DLAMCH 112* .. 113* .. Intrinsic Functions .. 114 INTRINSIC ABS, DBLE, MAX, MIN 115* .. 116* .. Executable Statements .. 117* 118 EPS = DLAMCH( 'Precision' ) 119 I1 = IFIRST 120 I2 = IFIRST 121 NEIG = ILAST - IFIRST + 1 122 NCNVRG = 0 123 THRESH = RELTOL 124 DO 10 I = IFIRST, ILAST 125 OPS = OPS + DBLE( 3 ) 126 IWORK( I ) = 0 127 PERT = EPS*( ABS( SIGMA )+ABS( W( I ) ) ) 128 WERR( I ) = WERR( I ) + PERT 129 IF( WGAP( I ).LT.PERT ) 130 $ WGAP( I ) = PERT 131 10 CONTINUE 132 DO 20 I = I1, ILAST 133 IF( I.EQ.1 ) THEN 134 GAP = WGAP( I ) 135 ELSE IF( I.EQ.N ) THEN 136 GAP = WGAP( I-1 ) 137 ELSE 138 GAP = MIN( WGAP( I-1 ), WGAP( I ) ) 139 END IF 140 OPS = OPS + DBLE( 1 ) 141 IF( WERR( I ).LT.THRESH*GAP ) THEN 142 NCNVRG = NCNVRG + 1 143 IWORK( I ) = 1 144 IF( I1.EQ.I ) 145 $ I1 = I1 + 1 146 ELSE 147 I2 = I 148 END IF 149 20 CONTINUE 150* 151* Initialize the unconverged intervals. 152* 153 I = I1 154 NINT = 0 155 RIGHT = ZERO 156 30 CONTINUE 157 IF( I.LE.I2 ) THEN 158 IF( IWORK( I ).EQ.0 ) THEN 159 DELTA = EPS 160 OPS = OPS + DBLE( 1 ) 161 LEFT = W( I ) - WERR( I ) 162* 163* Do while( CNT(LEFT).GT.I-1 ) 164* 165 40 CONTINUE 166 IF( I.GT.I1 .AND. LEFT.LE.RIGHT ) THEN 167 LEFT = RIGHT 168 CNT = I - 1 169 ELSE 170 S = -LEFT 171 CNT = 0 172 DO 50 J = 1, N - 1 173 OPS = OPS + DBLE( 5 ) 174 TMP = D( J ) + S 175 S = S*( LD( J ) / TMP )*L( J ) - LEFT 176 IF( TMP.LT.ZERO ) 177 $ CNT = CNT + 1 178 50 CONTINUE 179 TMP = D( N ) + S 180 IF( TMP.LT.ZERO ) 181 $ CNT = CNT + 1 182 IF( CNT.GT.I-1 ) THEN 183 OPS = OPS + DBLE( 4 ) 184 DELTA = TWO*DELTA 185 LEFT = LEFT - ( ABS( SIGMA )+ABS( LEFT ) )*DELTA 186 GO TO 40 187 END IF 188 END IF 189 OPS = OPS + DBLE( 1 ) 190 DELTA = EPS 191 RIGHT = W( I ) + WERR( I ) 192* 193* Do while( CNT(RIGHT).LT.I ) 194* 195 60 CONTINUE 196 S = -RIGHT 197 CNT = 0 198 OPS = OPS + DBLE( 5*( N-1 )+1 ) 199 DO 70 J = 1, N - 1 200 TMP = D( J ) + S 201 S = S*( LD( J ) / TMP )*L( J ) - RIGHT 202 IF( TMP.LT.ZERO ) 203 $ CNT = CNT + 1 204 70 CONTINUE 205 TMP = D( N ) + S 206 IF( TMP.LT.ZERO ) 207 $ CNT = CNT + 1 208 IF( CNT.LT.I ) THEN 209 OPS = OPS + DBLE( 4 ) 210 DELTA = TWO*DELTA 211 RIGHT = RIGHT + ( ABS( SIGMA )+ABS( RIGHT ) )*DELTA 212 GO TO 60 213 END IF 214 WERR( I ) = LEFT 215 W( I ) = RIGHT 216 IWORK( N+I ) = CNT 217 NINT = NINT + 1 218 I = CNT + 1 219 ELSE 220 I = I + 1 221 END IF 222 GO TO 30 223 END IF 224* 225* While( NCNVRG.LT.NEIG ) 226* 227 INITI1 = I1 228 INITI2 = I2 229 80 CONTINUE 230 IF( NCNVRG.LT.NEIG ) THEN 231 OLNINT = NINT 232 I = I1 233 DO 100 K = 1, OLNINT 234 NRIGHT = IWORK( N+I ) 235 IF( IWORK( I ).EQ.0 ) THEN 236 OPS = OPS + DBLE( 2 ) 237 MID = HALF*( WERR( I )+W( I ) ) 238 S = -MID 239 CNT = 0 240 OPS = OPS + DBLE( 5*( N-1 )+1 ) 241 DO 90 J = 1, N - 1 242 TMP = D( J ) + S 243 S = S*( LD( J ) / TMP )*L( J ) - MID 244 IF( TMP.LT.ZERO ) 245 $ CNT = CNT + 1 246 90 CONTINUE 247 TMP = D( N ) + S 248 IF( TMP.LT.ZERO ) 249 $ CNT = CNT + 1 250 CNT = MAX( I-1, MIN( NRIGHT, CNT ) ) 251 IF( I.EQ.NRIGHT ) THEN 252 IF( I.EQ.IFIRST ) THEN 253 OPS = OPS + DBLE( 1 ) 254 GAP = WERR( I+1 ) - W( I ) 255 ELSE IF( I.EQ.ILAST ) THEN 256 OPS = OPS + DBLE( 1 ) 257 GAP = WERR( I ) - W( I-1 ) 258 ELSE 259 OPS = OPS + DBLE( 2 ) 260 GAP = MIN( WERR( I+1 )-W( I ), WERR( I )-W( I-1 ) ) 261 END IF 262 OPS = OPS + DBLE( 2 ) 263 WIDTH = W( I ) - MID 264 IF( WIDTH.LT.THRESH*GAP ) THEN 265 NCNVRG = NCNVRG + 1 266 IWORK( I ) = 1 267 IF( I1.EQ.I ) THEN 268 I1 = I1 + 1 269 NINT = NINT - 1 270 END IF 271 END IF 272 END IF 273 IF( IWORK( I ).EQ.0 ) 274 $ I2 = K 275 IF( CNT.EQ.I-1 ) THEN 276 WERR( I ) = MID 277 ELSE IF( CNT.EQ.NRIGHT ) THEN 278 W( I ) = MID 279 ELSE 280 IWORK( N+I ) = CNT 281 NINT = NINT + 1 282 WERR( CNT+1 ) = MID 283 W( CNT+1 ) = W( I ) 284 W( I ) = MID 285 I = CNT + 1 286 IWORK( N+I ) = NRIGHT 287 END IF 288 END IF 289 I = NRIGHT + 1 290 100 CONTINUE 291 NINT = NINT - OLNINT + I2 292 GO TO 80 293 END IF 294 DO 110 I = INITI1, INITI2 295 OPS = OPS + DBLE( 3 ) 296 W( I ) = HALF*( WERR( I )+W( I ) ) 297 WERR( I ) = W( I ) - WERR( I ) 298 110 CONTINUE 299* 300 RETURN 301* 302* End of DLARRB 303* 304 END 305