github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dlasqtest/dlasq3.f (about) 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] Z 62 *> \verbatim 63 *> Z is DOUBLE PRECISION array, dimension ( 4*N ) 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[out] NFAIL 101 *> \verbatim 102 *> NFAIL is INTEGER 103 *> Number of times shift was too big. 104 *> \endverbatim 105 *> 106 *> \param[out] ITER 107 *> \verbatim 108 *> ITER is INTEGER 109 *> Number of iterations. 110 *> \endverbatim 111 *> 112 *> \param[out] NDIV 113 *> \verbatim 114 *> NDIV is INTEGER 115 *> Number of divisions. 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 September 2012 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.4.2) -- 186 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 187 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 188 * September 2012 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 227 N0IN = N0 228 EPS = DLAMCH( 'Precision' ) 229 TOL = EPS*HUNDRD 230 TOL2 = TOL**2 231 * 232 * Check for deflation. 233 * 234 10 CONTINUE 235 * 236 IF( N0.LT.I0 ) 237 $ RETURN 238 IF( N0.EQ.I0 ) 239 $ GO TO 20 240 NN = 4*N0 + PP 241 IF( N0.EQ.( I0+1 ) ) 242 $ GO TO 40 243 * 244 * Check whether E(N0-1) is negligible, 1 eigenvalue. 245 * 246 IF( Z( NN-5 ).GT.TOL2*( SIGMA+Z( NN-3 ) ) .AND. 247 $ Z( NN-2*PP-4 ).GT.TOL2*Z( NN-7 ) ) 248 $ GO TO 30 249 * 250 20 CONTINUE 251 * 252 Z( 4*N0-3 ) = Z( 4*N0+PP-3 ) + SIGMA 253 N0 = N0 - 1 254 GO TO 10 255 * 256 * Check whether E(N0-2) is negligible, 2 eigenvalues. 257 * 258 30 CONTINUE 259 * 260 IF( Z( NN-9 ).GT.TOL2*SIGMA .AND. 261 $ Z( NN-2*PP-8 ).GT.TOL2*Z( NN-11 ) ) 262 $ GO TO 50 263 * 264 40 CONTINUE 265 * 266 IF( Z( NN-3 ).GT.Z( NN-7 ) ) THEN 267 S = Z( NN-3 ) 268 Z( NN-3 ) = Z( NN-7 ) 269 Z( NN-7 ) = S 270 END IF 271 T = HALF*( ( Z( NN-7 )-Z( NN-3 ) )+Z( NN-5 ) ) 272 IF( Z( NN-5 ).GT.Z( NN-3 )*TOL2.AND.T.NE.ZERO ) THEN 273 S = Z( NN-3 )*( Z( NN-5 ) / T ) 274 IF( S.LE.T ) THEN 275 S = Z( NN-3 )*( Z( NN-5 ) / 276 $ ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 277 ELSE 278 S = Z( NN-3 )*( Z( NN-5 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 279 END IF 280 T = Z( NN-7 ) + ( S+Z( NN-5 ) ) 281 Z( NN-3 ) = Z( NN-3 )*( Z( NN-7 ) / T ) 282 Z( NN-7 ) = T 283 END IF 284 Z( 4*N0-7 ) = Z( NN-7 ) + SIGMA 285 Z( 4*N0-3 ) = Z( NN-3 ) + SIGMA 286 N0 = N0 - 2 287 GO TO 10 288 * 289 50 CONTINUE 290 IF( PP.EQ.2 ) 291 $ PP = 0 292 * 293 * Reverse the qd-array, if warranted. 294 * 295 296 IF( DMIN.LE.ZERO .OR. N0.LT.N0IN ) THEN 297 IF( CBIAS*Z( 4*I0+PP-3 ).LT.Z( 4*N0+PP-3 ) ) THEN 298 IPN4 = 4*( I0+N0 ) 299 DO 60 J4 = 4*I0, 2*( I0+N0-1 ), 4 300 TEMP = Z( J4-3 ) 301 Z( J4-3 ) = Z( IPN4-J4-3 ) 302 Z( IPN4-J4-3 ) = TEMP 303 TEMP = Z( J4-2 ) 304 Z( J4-2 ) = Z( IPN4-J4-2 ) 305 Z( IPN4-J4-2 ) = TEMP 306 TEMP = Z( J4-1 ) 307 Z( J4-1 ) = Z( IPN4-J4-5 ) 308 Z( IPN4-J4-5 ) = TEMP 309 TEMP = Z( J4 ) 310 Z( J4 ) = Z( IPN4-J4-4 ) 311 Z( IPN4-J4-4 ) = TEMP 312 60 CONTINUE 313 IF( N0-I0.LE.4 ) THEN 314 Z( 4*N0+PP-1 ) = Z( 4*I0+PP-1 ) 315 Z( 4*N0-PP ) = Z( 4*I0-PP ) 316 END IF 317 DMIN2 = MIN( DMIN2, Z( 4*N0+PP-1 ) ) 318 Z( 4*N0+PP-1 ) = MIN( Z( 4*N0+PP-1 ), Z( 4*I0+PP-1 ), 319 $ Z( 4*I0+PP+3 ) ) 320 Z( 4*N0-PP ) = MIN( Z( 4*N0-PP ), Z( 4*I0-PP ), 321 $ Z( 4*I0-PP+4 ) ) 322 QMAX = MAX( QMAX, Z( 4*I0+PP-3 ), Z( 4*I0+PP+1 ) ) 323 DMIN = -ZERO 324 END IF 325 END IF 326 * 327 * Choose a shift. 328 * 329 ! Print out DLASQ4 test cases 330 write(4,*) "{" 331 write(4,'(9999(g0))',advance="no") "z: []float64{" 332 do i = 1, NN 333 write (4,'(99999(e24.16,a))',advance="no") z(i), "," 334 end do 335 write (4,*) "}," 336 write (4,*) "i0: ", I0, "," 337 write (4,*) "n0: ", N0, "," 338 write (4,*) "pp: ", PP, "," 339 write (4,*) "n0in: ", N0IN, "," 340 write (4,*) "dmin: ", DMIN, "," 341 write (4,*) "dmin1:", DMIN1, "," 342 write (4,*) "dmin2:", DMIN2, "," 343 write (4,*) "dn: ", DN, "," 344 write (4,*) "dn1: ", DN1, "," 345 write (4,*) "dn2: ", DN2, "," 346 write (4,*) "tau: ", TAU, "," 347 write (4,*) "ttype: ", TTYPE, "," 348 write (4,*) "g: ", G, "," 349 CALL DLASQ4( I0, N0, Z, PP, N0IN, DMIN, DMIN1, DMIN2, DN, DN1, 350 $ DN2, TAU, TTYPE, G ) 351 352 write(4,'(9999(g0))',advance="no") "zOut: []float64{" 353 do i = 1, NN 354 write (4,'(99999(e24.16,a))',advance="no") z(i), "," 355 end do 356 write (4,*) "}," 357 write (4,*) "tauOut: ", TAU, "," 358 write (4,*) "ttypeOut: ", TTYPE, "," 359 write (4,*) "gOut: ", G, "," 360 write (4,*) "}," 361 362 * 363 * Call dqds until DMIN > 0. 364 * 365 70 CONTINUE 366 * 367 368 write(5,*) "{" 369 write(5,'(9999(g0))',advance="no") "z: []float64{" 370 do i = 1, NN 371 write (5,'(99999(e24.16,a))',advance="no") z(i), "," 372 end do 373 write (5,*) "}," 374 write (5,*) "i0: ", I0, "," 375 write (5,*) "n0: ", N0, "," 376 write (5,*) "pp: ", PP, "," 377 write (5,*) "tau: ", TAU, "," 378 write (5,*) "sigma: ", SIGMA, "," 379 write (5,*) "dmin: ", DMIN, "," 380 write (5,*) "dmin1:", DMIN1, "," 381 write (5,*) "dmin2:", DMIN2, "," 382 write (5,*) "dn: ", DN, "," 383 write (5,*) "dnm1: ", DN1, "," 384 write (5,*) "dnm2: ", DN2, "," 385 386 387 CALL DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, 388 $ DN1, DN2, IEEE, EPS ) 389 390 391 392 write (5,*) "i0Out: ", I0, "," 393 write (5,*) "n0Out: ", N0, "," 394 write (5,*) "ppOut: ", PP, "," 395 write (5,*) "tauOut: ", TAU, "," 396 write (5,*) "sigmaOut: ", SIGMA, "," 397 write (5,*) "dminOut: ", DMIN, "," 398 write (5,*) "dmin1Out:", DMIN1, "," 399 write (5,*) "dmin2Out:", DMIN2, "," 400 write (5,*) "dnOut: ", DN, "," 401 write (5,*) "dnm1Out: ", DN1, "," 402 write (5,*) "dnm2Out: ", DN2, "," 403 write (5,*) "}," 404 405 * 406 NDIV = NDIV + ( N0-I0+2 ) 407 408 ITER = ITER + 1 409 * 410 * Check status. 411 * 412 413 IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN 414 * 415 * Success. 416 * 417 GO TO 90 418 * 419 ELSE IF( DMIN.LT.ZERO .AND. DMIN1.GT.ZERO .AND. 420 $ Z( 4*( N0-1 )-PP ).LT.TOL*( SIGMA+DN1 ) .AND. 421 $ ABS( DN ).LT.TOL*SIGMA ) THEN 422 423 * 424 * Convergence hidden by negative DN. 425 * 426 Z( 4*( N0-1 )-PP+2 ) = ZERO 427 DMIN = ZERO 428 GO TO 90 429 ELSE IF( DMIN.LT.ZERO ) THEN 430 431 * 432 * TAU too big. Select new TAU and try again. 433 * 434 NFAIL = NFAIL + 1 435 IF( TTYPE.LT.-22 ) THEN 436 * 437 * Failed twice. Play it safe. 438 * 439 TAU = ZERO 440 ELSE IF( DMIN1.GT.ZERO ) THEN 441 * 442 * Late failure. Gives excellent shift. 443 * 444 TAU = ( TAU+DMIN )*( ONE-TWO*EPS ) 445 TTYPE = TTYPE - 11 446 ELSE 447 * 448 * Early failure. Divide by 4. 449 * 450 TAU = QURTR*TAU 451 TTYPE = TTYPE - 12 452 END IF 453 GO TO 70 454 ELSE IF( DISNAN( DMIN ) ) THEN 455 * 456 * NaN. 457 * 458 IF( TAU.EQ.ZERO ) THEN 459 GO TO 80 460 ELSE 461 TAU = ZERO 462 GO TO 70 463 END IF 464 ELSE 465 * 466 * Possible underflow. Play it safe. 467 * 468 GO TO 80 469 END IF 470 * 471 * Risk of underflow. 472 * 473 80 CONTINUE 474 475 CALL DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DN1, DN2 ) 476 477 478 NDIV = NDIV + ( N0-I0+2 ) 479 ITER = ITER + 1 480 TAU = ZERO 481 * 482 90 CONTINUE 483 484 IF( TAU.LT.SIGMA ) THEN 485 DESIG = DESIG + TAU 486 T = SIGMA + DESIG 487 DESIG = DESIG - ( T-SIGMA ) 488 ELSE 489 T = SIGMA + TAU 490 DESIG = SIGMA - ( T-TAU ) + DESIG 491 END IF 492 SIGMA = T 493 * 494 RETURN 495 * 496 * End of DLASQ3 497 * 498 END