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*> \date September 2012 140* 141*> \ingroup auxOTHERcomputational 142* 143* ===================================================================== 144 SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, 145 $ DN, DNM1, DNM2, IEEE, EPS ) 146* 147* -- LAPACK computational routine (version 3.4.2) -- 148* -- LAPACK is a software package provided by Univ. of Tennessee, -- 149* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 150* September 2012 151* 152* .. Scalar Arguments .. 153 LOGICAL IEEE 154 INTEGER I0, N0, PP 155 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, 156 $ SIGMA, EPS 157* .. 158* .. Array Arguments .. 159 DOUBLE PRECISION Z( * ) 160* .. 161* 162* ===================================================================== 163* 164* .. Parameter .. 165 DOUBLE PRECISION ZERO, HALF 166 PARAMETER ( ZERO = 0.0D0, HALF = 0.5 ) 167* .. 168* .. Local Scalars .. 169 INTEGER J4, J4P2 170 DOUBLE PRECISION D, EMIN, TEMP, DTHRESH 171* .. 172* .. Intrinsic Functions .. 173 INTRINSIC MIN 174* .. 175* .. Executable Statements .. 176* 177 IF( ( N0-I0-1 ).LE.0 ) 178 $ RETURN 179* 180 DTHRESH = EPS*(SIGMA+TAU) 181 IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO 182 IF( TAU.NE.ZERO ) THEN 183 J4 = 4*I0 + PP - 3 184 EMIN = Z( J4+4 ) 185 D = Z( J4 ) - TAU 186 DMIN = D 187 DMIN1 = -Z( J4 ) 188* 189 IF( IEEE ) THEN 190* 191* Code for IEEE arithmetic. 192* 193 IF( PP.EQ.0 ) THEN 194 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 195 Z( J4-2 ) = D + Z( J4-1 ) 196 TEMP = Z( J4+1 ) / Z( J4-2 ) 197 D = D*TEMP - TAU 198 DMIN = MIN( DMIN, D ) 199 Z( J4 ) = Z( J4-1 )*TEMP 200 EMIN = MIN( Z( J4 ), EMIN ) 201 10 CONTINUE 202 ELSE 203 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 204 Z( J4-3 ) = D + Z( J4 ) 205 TEMP = Z( J4+2 ) / Z( J4-3 ) 206 D = D*TEMP - TAU 207 DMIN = MIN( DMIN, D ) 208 Z( J4-1 ) = Z( J4 )*TEMP 209 EMIN = MIN( Z( J4-1 ), EMIN ) 210 20 CONTINUE 211 END IF 212* 213* Unroll last two steps. 214* 215 DNM2 = D 216 DMIN2 = DMIN 217 J4 = 4*( N0-2 ) - PP 218 J4P2 = J4 + 2*PP - 1 219 Z( J4-2 ) = DNM2 + Z( J4P2 ) 220 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 221 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 222 DMIN = MIN( DMIN, DNM1 ) 223* 224 DMIN1 = DMIN 225 J4 = J4 + 4 226 J4P2 = J4 + 2*PP - 1 227 Z( J4-2 ) = DNM1 + Z( J4P2 ) 228 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 229 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 230 DMIN = MIN( DMIN, DN ) 231* 232 ELSE 233* 234* Code for non IEEE arithmetic. 235* 236 IF( PP.EQ.0 ) THEN 237 DO 30 J4 = 4*I0, 4*( N0-3 ), 4 238 Z( J4-2 ) = D + Z( J4-1 ) 239 IF( D.LT.ZERO ) THEN 240 RETURN 241 ELSE 242 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 243 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 244 END IF 245 DMIN = MIN( DMIN, D ) 246 EMIN = MIN( EMIN, Z( J4 ) ) 247 30 CONTINUE 248 ELSE 249 DO 40 J4 = 4*I0, 4*( N0-3 ), 4 250 Z( J4-3 ) = D + Z( J4 ) 251 IF( D.LT.ZERO ) THEN 252 RETURN 253 ELSE 254 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 255 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 256 END IF 257 DMIN = MIN( DMIN, D ) 258 EMIN = MIN( EMIN, Z( J4-1 ) ) 259 40 CONTINUE 260 END IF 261* 262* Unroll last two steps. 263* 264 DNM2 = D 265 DMIN2 = DMIN 266 J4 = 4*( N0-2 ) - PP 267 J4P2 = J4 + 2*PP - 1 268 Z( J4-2 ) = DNM2 + Z( J4P2 ) 269 IF( DNM2.LT.ZERO ) THEN 270 RETURN 271 ELSE 272 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 273 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 274 END IF 275 DMIN = MIN( DMIN, DNM1 ) 276* 277 DMIN1 = DMIN 278 J4 = J4 + 4 279 J4P2 = J4 + 2*PP - 1 280 Z( J4-2 ) = DNM1 + Z( J4P2 ) 281 IF( DNM1.LT.ZERO ) THEN 282 RETURN 283 ELSE 284 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 285 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 286 END IF 287 DMIN = MIN( DMIN, DN ) 288* 289 END IF 290 ELSE 291* This is the version that sets d's to zero if they are small enough 292 J4 = 4*I0 + PP - 3 293 EMIN = Z( J4+4 ) 294 D = Z( J4 ) - TAU 295 DMIN = D 296 DMIN1 = -Z( J4 ) 297 IF( IEEE ) THEN 298* 299* Code for IEEE arithmetic. 300* 301 IF( PP.EQ.0 ) THEN 302 DO 50 J4 = 4*I0, 4*( N0-3 ), 4 303 Z( J4-2 ) = D + Z( J4-1 ) 304 TEMP = Z( J4+1 ) / Z( J4-2 ) 305 D = D*TEMP - TAU 306 IF( D.LT.DTHRESH ) D = ZERO 307 DMIN = MIN( DMIN, D ) 308 Z( J4 ) = Z( J4-1 )*TEMP 309 EMIN = MIN( Z( J4 ), EMIN ) 310 50 CONTINUE 311 ELSE 312 DO 60 J4 = 4*I0, 4*( N0-3 ), 4 313 Z( J4-3 ) = D + Z( J4 ) 314 TEMP = Z( J4+2 ) / Z( J4-3 ) 315 D = D*TEMP - TAU 316 IF( D.LT.DTHRESH ) D = ZERO 317 DMIN = MIN( DMIN, D ) 318 Z( J4-1 ) = Z( J4 )*TEMP 319 EMIN = MIN( Z( J4-1 ), EMIN ) 320 60 CONTINUE 321 END IF 322* 323* Unroll last two steps. 324* 325 DNM2 = D 326 DMIN2 = DMIN 327 J4 = 4*( N0-2 ) - PP 328 J4P2 = J4 + 2*PP - 1 329 Z( J4-2 ) = DNM2 + Z( J4P2 ) 330 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 331 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 332 DMIN = MIN( DMIN, DNM1 ) 333* 334 DMIN1 = DMIN 335 J4 = J4 + 4 336 J4P2 = J4 + 2*PP - 1 337 Z( J4-2 ) = DNM1 + Z( J4P2 ) 338 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 339 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 340 DMIN = MIN( DMIN, DN ) 341* 342 ELSE 343* 344* Code for non IEEE arithmetic. 345* 346 IF( PP.EQ.0 ) THEN 347 DO 70 J4 = 4*I0, 4*( N0-3 ), 4 348 Z( J4-2 ) = D + Z( J4-1 ) 349 IF( D.LT.ZERO ) THEN 350 RETURN 351 ELSE 352 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 353 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 354 END IF 355 IF( D.LT.DTHRESH) D = ZERO 356 DMIN = MIN( DMIN, D ) 357 EMIN = MIN( EMIN, Z( J4 ) ) 358 70 CONTINUE 359 ELSE 360 DO 80 J4 = 4*I0, 4*( N0-3 ), 4 361 Z( J4-3 ) = D + Z( J4 ) 362 IF( D.LT.ZERO ) THEN 363 RETURN 364 ELSE 365 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 366 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 367 END IF 368 IF( D.LT.DTHRESH) D = ZERO 369 DMIN = MIN( DMIN, D ) 370 EMIN = MIN( EMIN, Z( J4-1 ) ) 371 80 CONTINUE 372 END IF 373* 374* Unroll last two steps. 375* 376 DNM2 = D 377 DMIN2 = DMIN 378 J4 = 4*( N0-2 ) - PP 379 J4P2 = J4 + 2*PP - 1 380 Z( J4-2 ) = DNM2 + Z( J4P2 ) 381 IF( DNM2.LT.ZERO ) THEN 382 RETURN 383 ELSE 384 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 385 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 386 END IF 387 DMIN = MIN( DMIN, DNM1 ) 388* 389 DMIN1 = DMIN 390 J4 = J4 + 4 391 J4P2 = J4 + 2*PP - 1 392 Z( J4-2 ) = DNM1 + Z( J4P2 ) 393 IF( DNM1.LT.ZERO ) THEN 394 RETURN 395 ELSE 396 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 397 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 398 END IF 399 DMIN = MIN( DMIN, DN ) 400* 401 END IF 402 END IF 403* 404 Z( J4+2 ) = DN 405 Z( 4*N0-PP ) = EMIN 406 RETURN 407* 408* End of DLASQ5 409* 410 END 411c $Id$ 412