github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dlasqtest/dlasq5.f (about) 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 178 IF( ( N0-I0-1 ).LE.0 ) 179 $ RETURN 180 * 181 DTHRESH = EPS*(SIGMA+TAU) 182 IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO 183 IF( TAU.NE.ZERO ) THEN 184 J4 = 4*I0 + PP - 3 185 EMIN = Z( J4+4 ) 186 D = Z( J4 ) - TAU 187 DMIN = D 188 DMIN1 = -Z( J4 ) 189 * 190 IF( IEEE ) THEN 191 * 192 * Code for IEEE arithmetic. 193 * 194 IF( PP.EQ.0 ) THEN 195 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 196 Z( J4-2 ) = D + Z( J4-1 ) 197 TEMP = Z( J4+1 ) / Z( J4-2 ) 198 D = D*TEMP - TAU 199 DMIN = MIN( DMIN, D ) 200 Z( J4 ) = Z( J4-1 )*TEMP 201 EMIN = MIN( Z( J4 ), EMIN ) 202 10 CONTINUE 203 ELSE 204 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 205 Z( J4-3 ) = D + Z( J4 ) 206 TEMP = Z( J4+2 ) / Z( J4-3 ) 207 D = D*TEMP - TAU 208 DMIN = MIN( DMIN, D ) 209 Z( J4-1 ) = Z( J4 )*TEMP 210 EMIN = MIN( Z( J4-1 ), EMIN ) 211 20 CONTINUE 212 END IF 213 214 * 215 * Unroll last two steps. 216 * 217 DNM2 = D 218 DMIN2 = DMIN 219 J4 = 4*( N0-2 ) - PP 220 J4P2 = J4 + 2*PP - 1 221 Z( J4-2 ) = DNM2 + Z( J4P2 ) 222 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 223 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 224 DMIN = MIN( DMIN, DNM1 ) 225 * 226 DMIN1 = DMIN 227 J4 = J4 + 4 228 J4P2 = J4 + 2*PP - 1 229 Z( J4-2 ) = DNM1 + Z( J4P2 ) 230 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 231 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 232 DMIN = MIN( DMIN, DN ) 233 * 234 ELSE 235 * 236 * Code for non IEEE arithmetic. 237 * 238 IF( PP.EQ.0 ) THEN 239 DO 30 J4 = 4*I0, 4*( N0-3 ), 4 240 Z( J4-2 ) = D + Z( J4-1 ) 241 IF( D.LT.ZERO ) THEN 242 RETURN 243 ELSE 244 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 245 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 246 END IF 247 DMIN = MIN( DMIN, D ) 248 EMIN = MIN( EMIN, Z( J4 ) ) 249 30 CONTINUE 250 ELSE 251 DO 40 J4 = 4*I0, 4*( N0-3 ), 4 252 Z( J4-3 ) = D + Z( J4 ) 253 IF( D.LT.ZERO ) THEN 254 RETURN 255 ELSE 256 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 257 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 258 END IF 259 DMIN = MIN( DMIN, D ) 260 EMIN = MIN( EMIN, Z( J4-1 ) ) 261 40 CONTINUE 262 END IF 263 * 264 * Unroll last two steps. 265 * 266 DNM2 = D 267 DMIN2 = DMIN 268 J4 = 4*( N0-2 ) - PP 269 J4P2 = J4 + 2*PP - 1 270 Z( J4-2 ) = DNM2 + Z( J4P2 ) 271 IF( DNM2.LT.ZERO ) THEN 272 RETURN 273 ELSE 274 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 275 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 276 END IF 277 DMIN = MIN( DMIN, DNM1 ) 278 * 279 DMIN1 = DMIN 280 J4 = J4 + 4 281 J4P2 = J4 + 2*PP - 1 282 Z( J4-2 ) = DNM1 + Z( J4P2 ) 283 IF( DNM1.LT.ZERO ) THEN 284 RETURN 285 ELSE 286 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 287 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 288 END IF 289 DMIN = MIN( DMIN, DN ) 290 * 291 END IF 292 ELSE 293 * This is the version that sets d's to zero if they are small enough 294 J4 = 4*I0 + PP - 3 295 EMIN = Z( J4+4 ) 296 D = Z( J4 ) - TAU 297 DMIN = D 298 DMIN1 = -Z( J4 ) 299 IF( IEEE ) THEN 300 * 301 * Code for IEEE arithmetic. 302 * 303 304 IF( PP.EQ.0 ) THEN 305 DO 50 J4 = 4*I0, 4*( N0-3 ), 4 306 Z( J4-2 ) = D + Z( J4-1 ) 307 TEMP = Z( J4+1 ) / Z( J4-2 ) 308 D = D*TEMP - TAU 309 IF( D.LT.DTHRESH ) D = ZERO 310 DMIN = MIN( DMIN, D ) 311 Z( J4 ) = Z( J4-1 )*TEMP 312 EMIN = MIN( Z( J4 ), EMIN ) 313 50 CONTINUE 314 ELSE 315 DO 60 J4 = 4*I0, 4*( N0-3 ), 4 316 Z( J4-3 ) = D + Z( J4 ) 317 TEMP = Z( J4+2 ) / Z( J4-3 ) 318 D = D*TEMP - TAU 319 IF( D.LT.DTHRESH ) D = ZERO 320 DMIN = MIN( DMIN, D ) 321 Z( J4-1 ) = Z( J4 )*TEMP 322 EMIN = MIN( Z( J4-1 ), EMIN ) 323 60 CONTINUE 324 END IF 325 * 326 * Unroll last two steps. 327 * 328 DNM2 = D 329 DMIN2 = DMIN 330 J4 = 4*( N0-2 ) - PP 331 J4P2 = J4 + 2*PP - 1 332 Z( J4-2 ) = DNM2 + Z( J4P2 ) 333 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 334 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 335 DMIN = MIN( DMIN, DNM1 ) 336 * 337 DMIN1 = DMIN 338 J4 = J4 + 4 339 J4P2 = J4 + 2*PP - 1 340 Z( J4-2 ) = DNM1 + Z( J4P2 ) 341 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 342 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 343 DMIN = MIN( DMIN, DN ) 344 * 345 ELSE 346 * 347 * Code for non IEEE arithmetic. 348 * 349 IF( PP.EQ.0 ) THEN 350 DO 70 J4 = 4*I0, 4*( N0-3 ), 4 351 Z( J4-2 ) = D + Z( J4-1 ) 352 IF( D.LT.ZERO ) THEN 353 RETURN 354 ELSE 355 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 356 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 357 END IF 358 IF( D.LT.DTHRESH) D = ZERO 359 DMIN = MIN( DMIN, D ) 360 EMIN = MIN( EMIN, Z( J4 ) ) 361 70 CONTINUE 362 ELSE 363 DO 80 J4 = 4*I0, 4*( N0-3 ), 4 364 Z( J4-3 ) = D + Z( J4 ) 365 IF( D.LT.ZERO ) THEN 366 RETURN 367 ELSE 368 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 369 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 370 END IF 371 IF( D.LT.DTHRESH) D = ZERO 372 DMIN = MIN( DMIN, D ) 373 EMIN = MIN( EMIN, Z( J4-1 ) ) 374 80 CONTINUE 375 END IF 376 * 377 * Unroll last two steps. 378 * 379 DNM2 = D 380 DMIN2 = DMIN 381 J4 = 4*( N0-2 ) - PP 382 J4P2 = J4 + 2*PP - 1 383 Z( J4-2 ) = DNM2 + Z( J4P2 ) 384 IF( DNM2.LT.ZERO ) THEN 385 RETURN 386 ELSE 387 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 388 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 389 END IF 390 DMIN = MIN( DMIN, DNM1 ) 391 * 392 DMIN1 = DMIN 393 J4 = J4 + 4 394 J4P2 = J4 + 2*PP - 1 395 Z( J4-2 ) = DNM1 + Z( J4P2 ) 396 IF( DNM1.LT.ZERO ) THEN 397 RETURN 398 ELSE 399 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 400 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 401 END IF 402 DMIN = MIN( DMIN, DN ) 403 * 404 END IF 405 END IF 406 * 407 Z( J4+2 ) = DN 408 Z( 4*N0-PP ) = EMIN 409 RETURN 410 * 411 * End of DLASQ5 412 * 413 END