1*> \brief \b DLASQ4 computes an approximation to the smallest eigenvalue using values of d from the previous transform. Used by sbdsqr. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLASQ4 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq4.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq4.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq4.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, 22* DN1, DN2, TAU, TTYPE, G ) 23* 24* .. Scalar Arguments .. 25* INTEGER I0, N0, N0IN, PP, TTYPE 26* DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU 27* .. 28* .. Array Arguments .. 29* DOUBLE PRECISION Z( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> DLASQ4 computes an approximation TAU to the smallest eigenvalue 39*> using values of d from the previous transform. 40*> \endverbatim 41* 42* Arguments: 43* ========== 44* 45*> \param[in] I0 46*> \verbatim 47*> I0 is INTEGER 48*> First index. 49*> \endverbatim 50*> 51*> \param[in] N0 52*> \verbatim 53*> N0 is INTEGER 54*> Last index. 55*> \endverbatim 56*> 57*> \param[in] Z 58*> \verbatim 59*> Z is DOUBLE PRECISION array, dimension ( 4*N0 ) 60*> Z holds the qd array. 61*> \endverbatim 62*> 63*> \param[in] PP 64*> \verbatim 65*> PP is INTEGER 66*> PP=0 for ping, PP=1 for pong. 67*> \endverbatim 68*> 69*> \param[in] N0IN 70*> \verbatim 71*> N0IN is INTEGER 72*> The value of N0 at start of EIGTEST. 73*> \endverbatim 74*> 75*> \param[in] DMIN 76*> \verbatim 77*> DMIN is DOUBLE PRECISION 78*> Minimum value of d. 79*> \endverbatim 80*> 81*> \param[in] DMIN1 82*> \verbatim 83*> DMIN1 is DOUBLE PRECISION 84*> Minimum value of d, excluding D( N0 ). 85*> \endverbatim 86*> 87*> \param[in] DMIN2 88*> \verbatim 89*> DMIN2 is DOUBLE PRECISION 90*> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 91*> \endverbatim 92*> 93*> \param[in] DN 94*> \verbatim 95*> DN is DOUBLE PRECISION 96*> d(N) 97*> \endverbatim 98*> 99*> \param[in] DN1 100*> \verbatim 101*> DN1 is DOUBLE PRECISION 102*> d(N-1) 103*> \endverbatim 104*> 105*> \param[in] DN2 106*> \verbatim 107*> DN2 is DOUBLE PRECISION 108*> d(N-2) 109*> \endverbatim 110*> 111*> \param[out] TAU 112*> \verbatim 113*> TAU is DOUBLE PRECISION 114*> This is the shift. 115*> \endverbatim 116*> 117*> \param[out] TTYPE 118*> \verbatim 119*> TTYPE is INTEGER 120*> Shift type. 121*> \endverbatim 122*> 123*> \param[in,out] G 124*> \verbatim 125*> G is DOUBLE PRECISION 126*> G is passed as an argument in order to save its value between 127*> calls to DLASQ4. 128*> \endverbatim 129* 130* Authors: 131* ======== 132* 133*> \author Univ. of Tennessee 134*> \author Univ. of California Berkeley 135*> \author Univ. of Colorado Denver 136*> \author NAG Ltd. 137* 138*> \date June 2016 139* 140*> \ingroup auxOTHERcomputational 141* 142*> \par Further Details: 143* ===================== 144*> 145*> \verbatim 146*> 147*> CNST1 = 9/16 148*> \endverbatim 149*> 150* ===================================================================== 151 SUBROUTINE DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, 152 $ DN1, DN2, TAU, TTYPE, G ) 153* 154* -- LAPACK computational routine (version 3.7.1) -- 155* -- LAPACK is a software package provided by Univ. of Tennessee, -- 156* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 157* June 2016 158* 159* .. Scalar Arguments .. 160 INTEGER I0, N0, N0IN, PP, TTYPE 161 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, TAU 162* .. 163* .. Array Arguments .. 164 DOUBLE PRECISION Z( * ) 165* .. 166* 167* ===================================================================== 168* 169* .. Parameters .. 170 DOUBLE PRECISION CNST1, CNST2, CNST3 171 PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0, 172 $ CNST3 = 1.050D0 ) 173 DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD 174 PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0, 175 $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0, 176 $ TWO = 2.0D0, HUNDRD = 100.0D0 ) 177* .. 178* .. Local Scalars .. 179 INTEGER I4, NN, NP 180 DOUBLE PRECISION A2, B1, B2, GAM, GAP1, GAP2, S 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC MAX, MIN, SQRT 184* .. 185* .. Executable Statements .. 186* 187* A negative DMIN forces the shift to take that absolute value 188* TTYPE records the type of shift. 189* 190 IF( DMIN.LE.ZERO ) THEN 191 TAU = -DMIN 192 TTYPE = -1 193 RETURN 194 END IF 195* 196 NN = 4*N0 + PP 197 IF( N0IN.EQ.N0 ) THEN 198* 199* No eigenvalues deflated. 200* 201 IF( DMIN.EQ.DN .OR. DMIN.EQ.DN1 ) THEN 202* 203 B1 = SQRT( Z( NN-3 ) )*SQRT( Z( NN-5 ) ) 204 B2 = SQRT( Z( NN-7 ) )*SQRT( Z( NN-9 ) ) 205 A2 = Z( NN-7 ) + Z( NN-5 ) 206* 207* Cases 2 and 3. 208* 209 IF( DMIN.EQ.DN .AND. DMIN1.EQ.DN1 ) THEN 210 GAP2 = DMIN2 - A2 - DMIN2*QURTR 211 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2 ) THEN 212 GAP1 = A2 - DN - ( B2 / GAP2 )*B2 213 ELSE 214 GAP1 = A2 - DN - ( B1+B2 ) 215 END IF 216 IF( GAP1.GT.ZERO .AND. GAP1.GT.B1 ) THEN 217 S = MAX( DN-( B1 / GAP1 )*B1, HALF*DMIN ) 218 TTYPE = -2 219 ELSE 220 S = ZERO 221 IF( DN.GT.B1 ) 222 $ S = DN - B1 223 IF( A2.GT.( B1+B2 ) ) 224 $ S = MIN( S, A2-( B1+B2 ) ) 225 S = MAX( S, THIRD*DMIN ) 226 TTYPE = -3 227 END IF 228 ELSE 229* 230* Case 4. 231* 232 TTYPE = -4 233 S = QURTR*DMIN 234 IF( DMIN.EQ.DN ) THEN 235 GAM = DN 236 A2 = ZERO 237 IF( Z( NN-5 ) .GT. Z( NN-7 ) ) 238 $ RETURN 239 B2 = Z( NN-5 ) / Z( NN-7 ) 240 NP = NN - 9 241 ELSE 242 NP = NN - 2*PP 243 GAM = DN1 244 IF( Z( NP-4 ) .GT. Z( NP-2 ) ) 245 $ RETURN 246 A2 = Z( NP-4 ) / Z( NP-2 ) 247 IF( Z( NN-9 ) .GT. Z( NN-11 ) ) 248 $ RETURN 249 B2 = Z( NN-9 ) / Z( NN-11 ) 250 NP = NN - 13 251 END IF 252* 253* Approximate contribution to norm squared from I < NN-1. 254* 255 A2 = A2 + B2 256 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 257 IF( B2.EQ.ZERO ) 258 $ GO TO 20 259 B1 = B2 260 IF( Z( I4 ) .GT. Z( I4-2 ) ) 261 $ RETURN 262 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 263 A2 = A2 + B2 264 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 265 $ GO TO 20 266 10 CONTINUE 267 20 CONTINUE 268 A2 = CNST3*A2 269* 270* Rayleigh quotient residual bound. 271* 272 IF( A2.LT.CNST1 ) 273 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 274 END IF 275 ELSE IF( DMIN.EQ.DN2 ) THEN 276* 277* Case 5. 278* 279 TTYPE = -5 280 S = QURTR*DMIN 281* 282* Compute contribution to norm squared from I > NN-2. 283* 284 NP = NN - 2*PP 285 B1 = Z( NP-2 ) 286 B2 = Z( NP-6 ) 287 GAM = DN2 288 IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) 289 $ RETURN 290 A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) 291* 292* Approximate contribution to norm squared from I < NN-2. 293* 294 IF( N0-I0.GT.2 ) THEN 295 B2 = Z( NN-13 ) / Z( NN-15 ) 296 A2 = A2 + B2 297 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 298 IF( B2.EQ.ZERO ) 299 $ GO TO 40 300 B1 = B2 301 IF( Z( I4 ) .GT. Z( I4-2 ) ) 302 $ RETURN 303 B2 = B2*( Z( I4 ) / Z( I4-2 ) ) 304 A2 = A2 + B2 305 IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) 306 $ GO TO 40 307 30 CONTINUE 308 40 CONTINUE 309 A2 = CNST3*A2 310 END IF 311* 312 IF( A2.LT.CNST1 ) 313 $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) 314 ELSE 315* 316* Case 6, no information to guide us. 317* 318 IF( TTYPE.EQ.-6 ) THEN 319 G = G + THIRD*( ONE-G ) 320 ELSE IF( TTYPE.EQ.-18 ) THEN 321 G = QURTR*THIRD 322 ELSE 323 G = QURTR 324 END IF 325 S = G*DMIN 326 TTYPE = -6 327 END IF 328* 329 ELSE IF( N0IN.EQ.( N0+1 ) ) THEN 330* 331* One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. 332* 333 IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN 334* 335* Cases 7 and 8. 336* 337 TTYPE = -7 338 S = THIRD*DMIN1 339 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 340 $ RETURN 341 B1 = Z( NN-5 ) / Z( NN-7 ) 342 B2 = B1 343 IF( B2.EQ.ZERO ) 344 $ GO TO 60 345 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 346 A2 = B1 347 IF( Z( I4 ).GT.Z( I4-2 ) ) 348 $ RETURN 349 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 350 B2 = B2 + B1 351 IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) 352 $ GO TO 60 353 50 CONTINUE 354 60 CONTINUE 355 B2 = SQRT( CNST3*B2 ) 356 A2 = DMIN1 / ( ONE+B2**2 ) 357 GAP2 = HALF*DMIN2 - A2 358 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 359 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 360 ELSE 361 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 362 TTYPE = -8 363 END IF 364 ELSE 365* 366* Case 9. 367* 368 S = QURTR*DMIN1 369 IF( DMIN1.EQ.DN1 ) 370 $ S = HALF*DMIN1 371 TTYPE = -9 372 END IF 373* 374 ELSE IF( N0IN.EQ.( N0+2 ) ) THEN 375* 376* Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. 377* 378* Cases 10 and 11. 379* 380 IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN 381 TTYPE = -10 382 S = THIRD*DMIN2 383 IF( Z( NN-5 ).GT.Z( NN-7 ) ) 384 $ RETURN 385 B1 = Z( NN-5 ) / Z( NN-7 ) 386 B2 = B1 387 IF( B2.EQ.ZERO ) 388 $ GO TO 80 389 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 390 IF( Z( I4 ).GT.Z( I4-2 ) ) 391 $ RETURN 392 B1 = B1*( Z( I4 ) / Z( I4-2 ) ) 393 B2 = B2 + B1 394 IF( HUNDRD*B1.LT.B2 ) 395 $ GO TO 80 396 70 CONTINUE 397 80 CONTINUE 398 B2 = SQRT( CNST3*B2 ) 399 A2 = DMIN2 / ( ONE+B2**2 ) 400 GAP2 = Z( NN-7 ) + Z( NN-9 ) - 401 $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 402 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN 403 S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) 404 ELSE 405 S = MAX( S, A2*( ONE-CNST2*B2 ) ) 406 END IF 407 ELSE 408 S = QURTR*DMIN2 409 TTYPE = -11 410 END IF 411 ELSE IF( N0IN.GT.( N0+2 ) ) THEN 412* 413* Case 12, more than two eigenvalues deflated. No information. 414* 415 S = ZERO 416 TTYPE = -12 417 END IF 418* 419 TAU = S 420 RETURN 421* 422* End of DLASQ4 423* 424 END 425