1*> \brief <b> SLASQ5 computes one dqds transform in ping-pong form. Used by sbdsqr and sstegr. </b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLASQ5 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq5.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq5.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq5.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLASQ5( 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* REAL EPS, DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, SIGMA, TAU 28* .. 29* .. Array Arguments .. 30* REAL Z( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLASQ5 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 REAL 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 REAL 74*> This is the shift. 75*> \endverbatim 76*> 77*> \param[in] SIGMA 78*> \verbatim 79*> SIGMA is REAL 80*> This is the accumulated shift up to this step. 81*> \endverbatim 82*> 83*> \param[out] DMIN 84*> \verbatim 85*> DMIN is REAL 86*> Minimum value of d. 87*> \endverbatim 88*> 89*> \param[out] DMIN1 90*> \verbatim 91*> DMIN1 is REAL 92*> Minimum value of d, excluding D( N0 ). 93*> \endverbatim 94*> 95*> \param[out] DMIN2 96*> \verbatim 97*> DMIN2 is REAL 98*> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 99*> \endverbatim 100*> 101*> \param[out] DN 102*> \verbatim 103*> DN is REAL 104*> d(N0), the last value of d. 105*> \endverbatim 106*> 107*> \param[out] DNM1 108*> \verbatim 109*> DNM1 is REAL 110*> d(N0-1). 111*> \endverbatim 112*> 113*> \param[out] DNM2 114*> \verbatim 115*> DNM2 is REAL 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 REAL 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 December 2016 140* 141*> \ingroup auxOTHERcomputational 142* 143* ===================================================================== 144 SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, 145 $ DN, DNM1, DNM2, IEEE, EPS ) 146* 147* -- LAPACK computational routine (version 3.7.0) -- 148* -- LAPACK is a software package provided by Univ. of Tennessee, -- 149* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 150* December 2016 151* 152* .. Scalar Arguments .. 153 LOGICAL IEEE 154 INTEGER I0, N0, PP 155 REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, 156 $ SIGMA, EPS 157* .. 158* .. Array Arguments .. 159 REAL Z( * ) 160* .. 161* 162* ===================================================================== 163* 164* .. Parameter .. 165 REAL ZERO, HALF 166 PARAMETER ( ZERO = 0.0E0, HALF = 0.5 ) 167* .. 168* .. Local Scalars .. 169 INTEGER J4, J4P2 170 REAL 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* 291 ELSE 292* This is the version that sets d's to zero if they are small enough 293 J4 = 4*I0 + PP - 3 294 EMIN = Z( J4+4 ) 295 D = Z( J4 ) - TAU 296 DMIN = D 297 DMIN1 = -Z( J4 ) 298 IF( IEEE ) THEN 299* 300* Code for IEEE arithmetic. 301* 302 IF( PP.EQ.0 ) THEN 303 DO 50 J4 = 4*I0, 4*( N0-3 ), 4 304 Z( J4-2 ) = D + Z( J4-1 ) 305 TEMP = Z( J4+1 ) / Z( J4-2 ) 306 D = D*TEMP - TAU 307 IF( D.LT.DTHRESH ) D = ZERO 308 DMIN = MIN( DMIN, D ) 309 Z( J4 ) = Z( J4-1 )*TEMP 310 EMIN = MIN( Z( J4 ), EMIN ) 311 50 CONTINUE 312 ELSE 313 DO 60 J4 = 4*I0, 4*( N0-3 ), 4 314 Z( J4-3 ) = D + Z( J4 ) 315 TEMP = Z( J4+2 ) / Z( J4-3 ) 316 D = D*TEMP - TAU 317 IF( D.LT.DTHRESH ) D = ZERO 318 DMIN = MIN( DMIN, D ) 319 Z( J4-1 ) = Z( J4 )*TEMP 320 EMIN = MIN( Z( J4-1 ), EMIN ) 321 60 CONTINUE 322 END IF 323* 324* Unroll last two steps. 325* 326 DNM2 = D 327 DMIN2 = DMIN 328 J4 = 4*( N0-2 ) - PP 329 J4P2 = J4 + 2*PP - 1 330 Z( J4-2 ) = DNM2 + Z( J4P2 ) 331 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 332 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 333 DMIN = MIN( DMIN, DNM1 ) 334* 335 DMIN1 = DMIN 336 J4 = J4 + 4 337 J4P2 = J4 + 2*PP - 1 338 Z( J4-2 ) = DNM1 + Z( J4P2 ) 339 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 340 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 341 DMIN = MIN( DMIN, DN ) 342* 343 ELSE 344* 345* Code for non IEEE arithmetic. 346* 347 IF( PP.EQ.0 ) THEN 348 DO 70 J4 = 4*I0, 4*( N0-3 ), 4 349 Z( J4-2 ) = D + Z( J4-1 ) 350 IF( D.LT.ZERO ) THEN 351 RETURN 352 ELSE 353 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 354 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 355 END IF 356 IF( D.LT.DTHRESH ) D = ZERO 357 DMIN = MIN( DMIN, D ) 358 EMIN = MIN( EMIN, Z( J4 ) ) 359 70 CONTINUE 360 ELSE 361 DO 80 J4 = 4*I0, 4*( N0-3 ), 4 362 Z( J4-3 ) = D + Z( J4 ) 363 IF( D.LT.ZERO ) THEN 364 RETURN 365 ELSE 366 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 367 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 368 END IF 369 IF( D.LT.DTHRESH ) D = ZERO 370 DMIN = MIN( DMIN, D ) 371 EMIN = MIN( EMIN, Z( J4-1 ) ) 372 80 CONTINUE 373 END IF 374* 375* Unroll last two steps. 376* 377 DNM2 = D 378 DMIN2 = DMIN 379 J4 = 4*( N0-2 ) - PP 380 J4P2 = J4 + 2*PP - 1 381 Z( J4-2 ) = DNM2 + Z( J4P2 ) 382 IF( DNM2.LT.ZERO ) THEN 383 RETURN 384 ELSE 385 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 386 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 387 END IF 388 DMIN = MIN( DMIN, DNM1 ) 389* 390 DMIN1 = DMIN 391 J4 = J4 + 4 392 J4P2 = J4 + 2*PP - 1 393 Z( J4-2 ) = DNM1 + Z( J4P2 ) 394 IF( DNM1.LT.ZERO ) THEN 395 RETURN 396 ELSE 397 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 398 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 399 END IF 400 DMIN = MIN( DMIN, DN ) 401* 402 END IF 403* 404 END IF 405 Z( J4+2 ) = DN 406 Z( 4*N0-PP ) = EMIN 407 RETURN 408* 409* End of SLASQ5 410* 411 END 412