1*> \brief \b DLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLASQ5 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, 22* DNM1, DNM2, IEEE, EPS ) 23* 24* .. Scalar Arguments .. 25* LOGICAL IEEE 26* INTEGER I0, N0, PP 27* DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS 28* .. 29* .. Array Arguments .. 30* DOUBLE PRECISION Z( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DLASQ5 computes one dqds transform in ping-pong form, one 40*> version for IEEE machines another for non IEEE machines. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] I0 47*> \verbatim 48*> I0 is INTEGER 49*> First index. 50*> \endverbatim 51*> 52*> \param[in] N0 53*> \verbatim 54*> N0 is INTEGER 55*> Last index. 56*> \endverbatim 57*> 58*> \param[in] Z 59*> \verbatim 60*> Z is DOUBLE PRECISION array, dimension ( 4*N ) 61*> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid 62*> an extra argument. 63*> \endverbatim 64*> 65*> \param[in] PP 66*> \verbatim 67*> PP is INTEGER 68*> PP=0 for ping, PP=1 for pong. 69*> \endverbatim 70*> 71*> \param[in] TAU 72*> \verbatim 73*> TAU is DOUBLE PRECISION 74*> This is the shift. 75*> \endverbatim 76*> 77*> \param[in] SIGMA 78*> \verbatim 79*> SIGMA is DOUBLE PRECISION 80*> This is the accumulated shift up to this step. 81*> \endverbatim 82*> 83*> \param[out] DMIN 84*> \verbatim 85*> DMIN is DOUBLE PRECISION 86*> Minimum value of d. 87*> \endverbatim 88*> 89*> \param[out] DMIN1 90*> \verbatim 91*> DMIN1 is DOUBLE PRECISION 92*> Minimum value of d, excluding D( N0 ). 93*> \endverbatim 94*> 95*> \param[out] DMIN2 96*> \verbatim 97*> DMIN2 is DOUBLE PRECISION 98*> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 99*> \endverbatim 100*> 101*> \param[out] DN 102*> \verbatim 103*> DN is DOUBLE PRECISION 104*> d(N0), the last value of d. 105*> \endverbatim 106*> 107*> \param[out] DNM1 108*> \verbatim 109*> DNM1 is DOUBLE PRECISION 110*> d(N0-1). 111*> \endverbatim 112*> 113*> \param[out] DNM2 114*> \verbatim 115*> DNM2 is DOUBLE PRECISION 116*> d(N0-2). 117*> \endverbatim 118*> 119*> \param[in] IEEE 120*> \verbatim 121*> IEEE is LOGICAL 122*> Flag for IEEE or non IEEE arithmetic. 123*> \endverbatim 124*> 125*> \param[in] EPS 126*> \verbatim 127*> EPS is DOUBLE PRECISION 128*> This is the value of epsilon used. 129*> \endverbatim 130*> 131* Authors: 132* ======== 133* 134*> \author Univ. of Tennessee 135*> \author Univ. of California Berkeley 136*> \author Univ. of Colorado Denver 137*> \author NAG Ltd. 138* 139*> \ingroup auxOTHERcomputational 140* 141* ===================================================================== 142 SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, 143 $ DN, DNM1, DNM2, IEEE, EPS ) 144* 145* -- LAPACK computational routine -- 146* -- LAPACK is a software package provided by Univ. of Tennessee, -- 147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 148* 149* .. Scalar Arguments .. 150 LOGICAL IEEE 151 INTEGER I0, N0, PP 152 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, 153 $ SIGMA, EPS 154* .. 155* .. Array Arguments .. 156 DOUBLE PRECISION Z( * ) 157* .. 158* 159* ===================================================================== 160* 161* .. Parameter .. 162 DOUBLE PRECISION ZERO, HALF 163 PARAMETER ( ZERO = 0.0D0, HALF = 0.5 ) 164* .. 165* .. Local Scalars .. 166 INTEGER J4, J4P2 167 DOUBLE PRECISION D, EMIN, TEMP, DTHRESH 168* .. 169* .. Intrinsic Functions .. 170 INTRINSIC MIN 171* .. 172* .. Executable Statements .. 173* 174 IF( ( N0-I0-1 ).LE.0 ) 175 $ RETURN 176* 177 DTHRESH = EPS*(SIGMA+TAU) 178 IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO 179 IF( TAU.NE.ZERO ) THEN 180 J4 = 4*I0 + PP - 3 181 EMIN = Z( J4+4 ) 182 D = Z( J4 ) - TAU 183 DMIN = D 184 DMIN1 = -Z( J4 ) 185* 186 IF( IEEE ) THEN 187* 188* Code for IEEE arithmetic. 189* 190 IF( PP.EQ.0 ) THEN 191 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 192 Z( J4-2 ) = D + Z( J4-1 ) 193 TEMP = Z( J4+1 ) / Z( J4-2 ) 194 D = D*TEMP - TAU 195 DMIN = MIN( DMIN, D ) 196 Z( J4 ) = Z( J4-1 )*TEMP 197 EMIN = MIN( Z( J4 ), EMIN ) 198 10 CONTINUE 199 ELSE 200 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 201 Z( J4-3 ) = D + Z( J4 ) 202 TEMP = Z( J4+2 ) / Z( J4-3 ) 203 D = D*TEMP - TAU 204 DMIN = MIN( DMIN, D ) 205 Z( J4-1 ) = Z( J4 )*TEMP 206 EMIN = MIN( Z( J4-1 ), EMIN ) 207 20 CONTINUE 208 END IF 209* 210* Unroll last two steps. 211* 212 DNM2 = D 213 DMIN2 = DMIN 214 J4 = 4*( N0-2 ) - PP 215 J4P2 = J4 + 2*PP - 1 216 Z( J4-2 ) = DNM2 + Z( J4P2 ) 217 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 218 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 219 DMIN = MIN( DMIN, DNM1 ) 220* 221 DMIN1 = DMIN 222 J4 = J4 + 4 223 J4P2 = J4 + 2*PP - 1 224 Z( J4-2 ) = DNM1 + Z( J4P2 ) 225 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 226 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 227 DMIN = MIN( DMIN, DN ) 228* 229 ELSE 230* 231* Code for non IEEE arithmetic. 232* 233 IF( PP.EQ.0 ) THEN 234 DO 30 J4 = 4*I0, 4*( N0-3 ), 4 235 Z( J4-2 ) = D + Z( J4-1 ) 236 IF( D.LT.ZERO ) THEN 237 RETURN 238 ELSE 239 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 240 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 241 END IF 242 DMIN = MIN( DMIN, D ) 243 EMIN = MIN( EMIN, Z( J4 ) ) 244 30 CONTINUE 245 ELSE 246 DO 40 J4 = 4*I0, 4*( N0-3 ), 4 247 Z( J4-3 ) = D + Z( J4 ) 248 IF( D.LT.ZERO ) THEN 249 RETURN 250 ELSE 251 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 252 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 253 END IF 254 DMIN = MIN( DMIN, D ) 255 EMIN = MIN( EMIN, Z( J4-1 ) ) 256 40 CONTINUE 257 END IF 258* 259* Unroll last two steps. 260* 261 DNM2 = D 262 DMIN2 = DMIN 263 J4 = 4*( N0-2 ) - PP 264 J4P2 = J4 + 2*PP - 1 265 Z( J4-2 ) = DNM2 + Z( J4P2 ) 266 IF( DNM2.LT.ZERO ) THEN 267 RETURN 268 ELSE 269 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 270 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 271 END IF 272 DMIN = MIN( DMIN, DNM1 ) 273* 274 DMIN1 = DMIN 275 J4 = J4 + 4 276 J4P2 = J4 + 2*PP - 1 277 Z( J4-2 ) = DNM1 + Z( J4P2 ) 278 IF( DNM1.LT.ZERO ) THEN 279 RETURN 280 ELSE 281 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 282 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 283 END IF 284 DMIN = MIN( DMIN, DN ) 285* 286 END IF 287 ELSE 288* This is the version that sets d's to zero if they are small enough 289 J4 = 4*I0 + PP - 3 290 EMIN = Z( J4+4 ) 291 D = Z( J4 ) - TAU 292 DMIN = D 293 DMIN1 = -Z( J4 ) 294 IF( IEEE ) THEN 295* 296* Code for IEEE arithmetic. 297* 298 IF( PP.EQ.0 ) THEN 299 DO 50 J4 = 4*I0, 4*( N0-3 ), 4 300 Z( J4-2 ) = D + Z( J4-1 ) 301 TEMP = Z( J4+1 ) / Z( J4-2 ) 302 D = D*TEMP - TAU 303 IF( D.LT.DTHRESH ) D = ZERO 304 DMIN = MIN( DMIN, D ) 305 Z( J4 ) = Z( J4-1 )*TEMP 306 EMIN = MIN( Z( J4 ), EMIN ) 307 50 CONTINUE 308 ELSE 309 DO 60 J4 = 4*I0, 4*( N0-3 ), 4 310 Z( J4-3 ) = D + Z( J4 ) 311 TEMP = Z( J4+2 ) / Z( J4-3 ) 312 D = D*TEMP - TAU 313 IF( D.LT.DTHRESH ) D = ZERO 314 DMIN = MIN( DMIN, D ) 315 Z( J4-1 ) = Z( J4 )*TEMP 316 EMIN = MIN( Z( J4-1 ), EMIN ) 317 60 CONTINUE 318 END IF 319* 320* Unroll last two steps. 321* 322 DNM2 = D 323 DMIN2 = DMIN 324 J4 = 4*( N0-2 ) - PP 325 J4P2 = J4 + 2*PP - 1 326 Z( J4-2 ) = DNM2 + Z( J4P2 ) 327 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 328 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 329 DMIN = MIN( DMIN, DNM1 ) 330* 331 DMIN1 = DMIN 332 J4 = J4 + 4 333 J4P2 = J4 + 2*PP - 1 334 Z( J4-2 ) = DNM1 + Z( J4P2 ) 335 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 336 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 337 DMIN = MIN( DMIN, DN ) 338* 339 ELSE 340* 341* Code for non IEEE arithmetic. 342* 343 IF( PP.EQ.0 ) THEN 344 DO 70 J4 = 4*I0, 4*( N0-3 ), 4 345 Z( J4-2 ) = D + Z( J4-1 ) 346 IF( D.LT.ZERO ) THEN 347 RETURN 348 ELSE 349 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 350 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 351 END IF 352 IF( D.LT.DTHRESH) D = ZERO 353 DMIN = MIN( DMIN, D ) 354 EMIN = MIN( EMIN, Z( J4 ) ) 355 70 CONTINUE 356 ELSE 357 DO 80 J4 = 4*I0, 4*( N0-3 ), 4 358 Z( J4-3 ) = D + Z( J4 ) 359 IF( D.LT.ZERO ) THEN 360 RETURN 361 ELSE 362 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 363 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 364 END IF 365 IF( D.LT.DTHRESH) D = ZERO 366 DMIN = MIN( DMIN, D ) 367 EMIN = MIN( EMIN, Z( J4-1 ) ) 368 80 CONTINUE 369 END IF 370* 371* Unroll last two steps. 372* 373 DNM2 = D 374 DMIN2 = DMIN 375 J4 = 4*( N0-2 ) - PP 376 J4P2 = J4 + 2*PP - 1 377 Z( J4-2 ) = DNM2 + Z( J4P2 ) 378 IF( DNM2.LT.ZERO ) THEN 379 RETURN 380 ELSE 381 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 382 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 383 END IF 384 DMIN = MIN( DMIN, DNM1 ) 385* 386 DMIN1 = DMIN 387 J4 = J4 + 4 388 J4P2 = J4 + 2*PP - 1 389 Z( J4-2 ) = DNM1 + Z( J4P2 ) 390 IF( DNM1.LT.ZERO ) THEN 391 RETURN 392 ELSE 393 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 394 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 395 END IF 396 DMIN = MIN( DMIN, DN ) 397* 398 END IF 399 END IF 400* 401 Z( J4+2 ) = DN 402 Z( 4*N0-PP ) = EMIN 403 RETURN 404* 405* End of DLASQ5 406* 407 END 408