1*> \brief \b DLASQ3 checks for deflation, computes a shift and calls dqds. 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 DLASQ3 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq3.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq3.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq3.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 22* ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 23* DN2, G, TAU ) 24* 25* .. Scalar Arguments .. 26* LOGICAL IEEE 27* INTEGER I0, ITER, N0, NDIV, NFAIL, PP 28* DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, 29* $ QMAX, SIGMA, TAU 30* .. 31* .. Array Arguments .. 32* DOUBLE PRECISION Z( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> DLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. 42*> In case of failure it changes shifts, and tries again until output 43*> is positive. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] I0 50*> \verbatim 51*> I0 is INTEGER 52*> First index. 53*> \endverbatim 54*> 55*> \param[in,out] N0 56*> \verbatim 57*> N0 is INTEGER 58*> Last index. 59*> \endverbatim 60*> 61*> \param[in,out] Z 62*> \verbatim 63*> Z is DOUBLE PRECISION array, dimension ( 4*N0 ) 64*> Z holds the qd array. 65*> \endverbatim 66*> 67*> \param[in,out] PP 68*> \verbatim 69*> PP is INTEGER 70*> PP=0 for ping, PP=1 for pong. 71*> PP=2 indicates that flipping was applied to the Z array 72*> and that the initial tests for deflation should not be 73*> performed. 74*> \endverbatim 75*> 76*> \param[out] DMIN 77*> \verbatim 78*> DMIN is DOUBLE PRECISION 79*> Minimum value of d. 80*> \endverbatim 81*> 82*> \param[out] SIGMA 83*> \verbatim 84*> SIGMA is DOUBLE PRECISION 85*> Sum of shifts used in current segment. 86*> \endverbatim 87*> 88*> \param[in,out] DESIG 89*> \verbatim 90*> DESIG is DOUBLE PRECISION 91*> Lower order part of SIGMA 92*> \endverbatim 93*> 94*> \param[in] QMAX 95*> \verbatim 96*> QMAX is DOUBLE PRECISION 97*> Maximum value of q. 98*> \endverbatim 99*> 100*> \param[in,out] NFAIL 101*> \verbatim 102*> NFAIL is INTEGER 103*> Increment NFAIL by 1 each time the shift was too big. 104*> \endverbatim 105*> 106*> \param[in,out] ITER 107*> \verbatim 108*> ITER is INTEGER 109*> Increment ITER by 1 for each iteration. 110*> \endverbatim 111*> 112*> \param[in,out] NDIV 113*> \verbatim 114*> NDIV is INTEGER 115*> Increment NDIV by 1 for each division. 116*> \endverbatim 117*> 118*> \param[in] IEEE 119*> \verbatim 120*> IEEE is LOGICAL 121*> Flag for IEEE or non IEEE arithmetic (passed to DLASQ5). 122*> \endverbatim 123*> 124*> \param[in,out] TTYPE 125*> \verbatim 126*> TTYPE is INTEGER 127*> Shift type. 128*> \endverbatim 129*> 130*> \param[in,out] DMIN1 131*> \verbatim 132*> DMIN1 is DOUBLE PRECISION 133*> \endverbatim 134*> 135*> \param[in,out] DMIN2 136*> \verbatim 137*> DMIN2 is DOUBLE PRECISION 138*> \endverbatim 139*> 140*> \param[in,out] DN 141*> \verbatim 142*> DN is DOUBLE PRECISION 143*> \endverbatim 144*> 145*> \param[in,out] DN1 146*> \verbatim 147*> DN1 is DOUBLE PRECISION 148*> \endverbatim 149*> 150*> \param[in,out] DN2 151*> \verbatim 152*> DN2 is DOUBLE PRECISION 153*> \endverbatim 154*> 155*> \param[in,out] G 156*> \verbatim 157*> G is DOUBLE PRECISION 158*> \endverbatim 159*> 160*> \param[in,out] TAU 161*> \verbatim 162*> TAU is DOUBLE PRECISION 163*> 164*> These are passed as arguments in order to save their values 165*> between calls to DLASQ3. 166*> \endverbatim 167* 168* Authors: 169* ======== 170* 171*> \author Univ. of Tennessee 172*> \author Univ. of California Berkeley 173*> \author Univ. of Colorado Denver 174*> \author NAG Ltd. 175* 176*> \date June 2016 177* 178*> \ingroup auxOTHERcomputational 179* 180* ===================================================================== 181 SUBROUTINE DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 182 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 183 $ DN2, G, TAU ) 184* 185* -- LAPACK computational routine (version 3.7.0) -- 186* -- LAPACK is a software package provided by Univ. of Tennessee, -- 187* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 188* June 2016 189* 190* .. Scalar Arguments .. 191 LOGICAL IEEE 192 INTEGER I0, ITER, N0, NDIV, NFAIL, PP 193 DOUBLE PRECISION DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, G, 194 $ QMAX, SIGMA, TAU 195* .. 196* .. Array Arguments .. 197 DOUBLE PRECISION Z( * ) 198* .. 199* 200* ===================================================================== 201* 202* .. Parameters .. 203 DOUBLE PRECISION CBIAS 204 PARAMETER ( CBIAS = 1.50D0 ) 205 DOUBLE PRECISION ZERO, QURTR, HALF, ONE, TWO, HUNDRD 206 PARAMETER ( ZERO = 0.0D0, QURTR = 0.250D0, HALF = 0.5D0, 207 $ ONE = 1.0D0, TWO = 2.0D0, HUNDRD = 100.0D0 ) 208* .. 209* .. Local Scalars .. 210 INTEGER IPN4, J4, N0IN, NN, TTYPE 211 DOUBLE PRECISION EPS, S, T, TEMP, TOL, TOL2 212* .. 213* .. External Subroutines .. 214 EXTERNAL DLASQ4, DLASQ5, DLASQ6 215* .. 216* .. External Function .. 217 DOUBLE PRECISION DLAMCH 218 LOGICAL DISNAN 219 EXTERNAL DISNAN, DLAMCH 220* .. 221* .. Intrinsic Functions .. 222 INTRINSIC ABS, MAX, MIN, SQRT 223* .. 224* .. Executable Statements .. 225* 226 N0IN = N0 227 EPS = DLAMCH( 'Precision' ) 228 TOL = EPS*HUNDRD 229 TOL2 = TOL**2 230* 231* Check for deflation. 232* 233 10 CONTINUE 234* 235 IF( N0.LT.I0 ) 236 $ RETURN 237 IF( N0.EQ.I0 ) 238 $ GO TO 20 239 NN = 4*N0 + PP 240 IF( N0.EQ.( I0+1 ) ) 241 $ GO TO 40 242* 243* Check whether E(N0-1) is negligible, 1 eigenvalue. 244* 245 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND. 246 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) ) 247 $ GO TO 30 248* 249 20 CONTINUE 250* 251 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA 252 N0 = N0 - 1 253 GO TO 10 254* 255* Check whether E(N0-2) is negligible, 2 eigenvalues. 256* 257 30 CONTINUE 258* 259 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND. 260 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) ) 261 $ GO TO 50 262* 263 40 CONTINUE 264* 265 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN 266 S = Z( NN-3 ) 267 Z( NN-3 ) = Z( NN-7 ) 268 Z( NN-7 ) = S 269 END IF 270 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) 271 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN 272 S = Z( NN-3 )*( Z( NN-5 ) / T ) 273 IF( S.LE.T ) THEN 274 S = Z( NN-3 )*( Z( NN-5 ) / 275 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 276 ELSE 277 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 278 END IF 279 T = Z( NN-7 ) + ( S+Z( NN-5 ) ) 280 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T ) 281 Z( NN-7 ) = T 282 END IF 283 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA 284 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA 285 N0 = N0 - 2 286 GO TO 10 287* 288 50 CONTINUE 289 IF( PP.EQ.2 ) 290 $ PP = 0 291* 292* Reverse the qd-array, if warranted. 293* 294 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN 295 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN 296 IPN4 = 4*( I0+N0 ) 297 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4 298 TEMP = Z( J4-3 ) 299 Z( J4-3 ) = Z( IPN4-J4-3 ) 300 Z( IPN4-J4-3 ) = TEMP 301 TEMP = Z( J4-2 ) 302 Z( J4-2 ) = Z( IPN4-J4-2 ) 303 Z( IPN4-J4-2 ) = TEMP 304 TEMP = Z( J4-1 ) 305 Z( J4-1 ) = Z( IPN4-J4-5 ) 306 Z( IPN4-J4-5 ) = TEMP 307 TEMP = Z( J4 ) 308 Z( J4 ) = Z( IPN4-J4-4 ) 309 Z( IPN4-J4-4 ) = TEMP 310 60 CONTINUE 311 IF( N0-I0.LE.4 ) THEN 312 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 ) 313 Z( 4*N0-PP ) = Z( 4*I0-PP ) 314 END IF 315 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) ) 316 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ), 317 $ Z( 4*I0+PP+3 ) ) 318 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ), 319 $ Z( 4*I0-PP+4 ) ) 320 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) ) 321 DMIN = -ZERO 322 END IF 323 END IF 324* 325* Choose a shift. 326* 327 CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, 328 $ DN2, TAU, TTYPE, G ) 329* 330* Call dqds until DMIN > 0. 331* 332 70 CONTINUE 333* 334 CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, 335 $ DN1, DN2, IEEE, EPS ) 336* 337 NDIV = NDIV + ( N0-I0+2 ) 338 ITER = ITER + 1 339* 340* Check status. 341* 342 IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN 343* 344* Success. 345* 346 GO TO 90 347* 348 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 349 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. 350 $ ABS( DN ).LT.TOL*SIGMA ) THEN 351* 352* Convergence hidden by negative DN. 353* 354 Z( 4*( N0-1 )-PP+2 ) = ZERO 355 DMIN = ZERO 356 GO TO 90 357 ELSE IF( DMIN.LT.ZERO ) THEN 358* 359* TAU too big. Select new TAU and try again. 360* 361 NFAIL = NFAIL + 1 362 IF( TTYPE.LT.-22 ) THEN 363* 364* Failed twice. Play it safe. 365* 366 TAU = ZERO 367 ELSE IF( DMIN1.GT.ZERO ) THEN 368* 369* Late failure. Gives excellent shift. 370* 371 TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) 372 TTYPE = TTYPE - 11 373 ELSE 374* 375* Early failure. Divide by 4. 376* 377 TAU = QURTR*TAU 378 TTYPE = TTYPE - 12 379 END IF 380 GO TO 70 381 ELSE IF( DISNAN( DMIN ) ) THEN 382* 383* NaN. 384* 385 IF( TAU.EQ.ZERO ) THEN 386 GO TO 80 387 ELSE 388 TAU = ZERO 389 GO TO 70 390 END IF 391 ELSE 392* 393* Possible underflow. Play it safe. 394* 395 GO TO 80 396 END IF 397* 398* Risk of underflow. 399* 400 80 CONTINUE 401 CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) 402 NDIV = NDIV + ( N0-I0+2 ) 403 ITER = ITER + 1 404 TAU = ZERO 405* 406 90 CONTINUE 407 IF( TAU.LT.SIGMA ) THEN 408 DESIG = DESIG + TAU 409 T = SIGMA + DESIG 410 DESIG = DESIG - ( T-SIGMA ) 411 ELSE 412 T = SIGMA + TAU 413 DESIG = SIGMA - ( T-TAU ) + DESIG 414 END IF 415 SIGMA = T 416* 417 RETURN 418* 419* End of DLASQ3 420* 421 END 422