1*> \brief \b DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. 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 DLASQ2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLASQ2( N, Z, INFO ) 22* 23* .. Scalar Arguments .. 24* INTEGER INFO, N 25* .. 26* .. Array Arguments .. 27* DOUBLE PRECISION Z( * ) 28* .. 29* 30* 31*> \par Purpose: 32* ============= 33*> 34*> \verbatim 35*> 36*> DLASQ2 computes all the eigenvalues of the symmetric positive 37*> definite tridiagonal matrix associated with the qd array Z to high 38*> relative accuracy are computed to high relative accuracy, in the 39*> absence of denormalization, underflow and overflow. 40*> 41*> To see the relation of Z to the tridiagonal matrix, let L be a 42*> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and 43*> let U be an upper bidiagonal matrix with 1's above and diagonal 44*> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the 45*> symmetric tridiagonal to which it is similar. 46*> 47*> Note : DLASQ2 defines a logical variable, IEEE, which is true 48*> on machines which follow ieee-754 floating-point standard in their 49*> handling of infinities and NaNs, and false otherwise. This variable 50*> is passed to DLASQ3. 51*> \endverbatim 52* 53* Arguments: 54* ========== 55* 56*> \param[in] N 57*> \verbatim 58*> N is INTEGER 59*> The number of rows and columns in the matrix. N >= 0. 60*> \endverbatim 61*> 62*> \param[in,out] Z 63*> \verbatim 64*> Z is DOUBLE PRECISION array, dimension ( 4*N ) 65*> On entry Z holds the qd array. On exit, entries 1 to N hold 66*> the eigenvalues in decreasing order, Z( 2*N+1 ) holds the 67*> trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If 68*> N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) 69*> holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of 70*> shifts that failed. 71*> \endverbatim 72*> 73*> \param[out] INFO 74*> \verbatim 75*> INFO is INTEGER 76*> = 0: successful exit 77*> < 0: if the i-th argument is a scalar and had an illegal 78*> value, then INFO = -i, if the i-th argument is an 79*> array and the j-entry had an illegal value, then 80*> INFO = -(i*100+j) 81*> > 0: the algorithm failed 82*> = 1, a split was marked by a positive value in E 83*> = 2, current block of Z not diagonalized after 100*N 84*> iterations (in inner while loop). On exit Z holds 85*> a qd array with the same eigenvalues as the given Z. 86*> = 3, termination criterion of outer while loop not met 87*> (program created more than N unreduced blocks) 88*> \endverbatim 89* 90* Authors: 91* ======== 92* 93*> \author Univ. of Tennessee 94*> \author Univ. of California Berkeley 95*> \author Univ. of Colorado Denver 96*> \author NAG Ltd. 97* 98*> \date December 2016 99* 100*> \ingroup auxOTHERcomputational 101* 102*> \par Further Details: 103* ===================== 104*> 105*> \verbatim 106*> 107*> Local Variables: I0:N0 defines a current unreduced segment of Z. 108*> The shifts are accumulated in SIGMA. Iteration count is in ITER. 109*> Ping-pong is controlled by PP (alternates between 0 and 1). 110*> \endverbatim 111*> 112* ===================================================================== 113 SUBROUTINE DLASQ2( N, Z, INFO ) 114* 115* -- LAPACK computational routine (version 3.7.0) -- 116* -- LAPACK is a software package provided by Univ. of Tennessee, -- 117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 118* December 2016 119* 120* .. Scalar Arguments .. 121 INTEGER INFO, N 122* .. 123* .. Array Arguments .. 124 DOUBLE PRECISION Z( * ) 125* .. 126* 127* ===================================================================== 128* 129* .. Parameters .. 130 DOUBLE PRECISION CBIAS 131 PARAMETER ( CBIAS = 1.50D0 ) 132 DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD 133 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, 134 $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 ) 135* .. 136* .. Local Scalars .. 137 LOGICAL IEEE 138 INTEGER I0, I1, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, 139 $ K, KMIN, N0, N1, NBIG, NDIV, NFAIL, PP, SPLT, 140 $ TTYPE 141 DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN, 142 $ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX, 143 $ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL, 144 $ TOL2, TRACE, ZMAX, TEMPE, TEMPQ 145* .. 146* .. External Subroutines .. 147 EXTERNAL DLASQ3, DLASRT, XERBLA 148* .. 149* .. External Functions .. 150 INTEGER ILAENV 151 DOUBLE PRECISION DLAMCH 152 EXTERNAL DLAMCH, ILAENV 153* .. 154* .. Intrinsic Functions .. 155 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 156* .. 157* .. Executable Statements .. 158* 159* Test the input arguments. 160* (in case DLASQ2 is not called by DLASQ1) 161* 162 INFO = 0 163 EPS = DLAMCH( 'Precision' ) 164 SAFMIN = DLAMCH( 'Safe minimum' ) 165 TOL = EPS*HUNDRD 166 TOL2 = TOL**2 167* 168 IF( N.LT.0 ) THEN 169 INFO = -1 170 CALL XERBLA( 'DLASQ2', 1 ) 171 RETURN 172 ELSE IF( N.EQ.0 ) THEN 173 RETURN 174 ELSE IF( N.EQ.1 ) THEN 175* 176* 1-by-1 case. 177* 178 IF( Z( 1 ).LT.ZERO ) THEN 179 INFO = -201 180 CALL XERBLA( 'DLASQ2', 2 ) 181 END IF 182 RETURN 183 ELSE IF( N.EQ.2 ) THEN 184* 185* 2-by-2 case. 186* 187 IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN 188 INFO = -2 189 CALL XERBLA( 'DLASQ2', 2 ) 190 RETURN 191 ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN 192 D = Z( 3 ) 193 Z( 3 ) = Z( 1 ) 194 Z( 1 ) = D 195 END IF 196 Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) 197 IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN 198 T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 199 S = Z( 3 )*( Z( 2 ) / T ) 200 IF( S.LE.T ) THEN 201 S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 202 ELSE 203 S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 204 END IF 205 T = Z( 1 ) + ( S+Z( 2 ) ) 206 Z( 3 ) = Z( 3 )*( Z( 1 ) / T ) 207 Z( 1 ) = T 208 END IF 209 Z( 2 ) = Z( 3 ) 210 Z( 6 ) = Z( 2 ) + Z( 1 ) 211 RETURN 212 END IF 213* 214* Check for negative data and compute sums of q's and e's. 215* 216 Z( 2*N ) = ZERO 217 EMIN = Z( 2 ) 218 QMAX = ZERO 219 ZMAX = ZERO 220 D = ZERO 221 E = ZERO 222* 223 DO 10 K = 1, 2*( N-1 ), 2 224 IF( Z( K ).LT.ZERO ) THEN 225 INFO = -( 200+K ) 226 CALL XERBLA( 'DLASQ2', 2 ) 227 RETURN 228 ELSE IF( Z( K+1 ).LT.ZERO ) THEN 229 INFO = -( 200+K+1 ) 230 CALL XERBLA( 'DLASQ2', 2 ) 231 RETURN 232 END IF 233 D = D + Z( K ) 234 E = E + Z( K+1 ) 235 QMAX = MAX( QMAX, Z( K ) ) 236 EMIN = MIN( EMIN, Z( K+1 ) ) 237 ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) ) 238 10 CONTINUE 239 IF( Z( 2*N-1 ).LT.ZERO ) THEN 240 INFO = -( 200+2*N-1 ) 241 CALL XERBLA( 'DLASQ2', 2 ) 242 RETURN 243 END IF 244 D = D + Z( 2*N-1 ) 245 QMAX = MAX( QMAX, Z( 2*N-1 ) ) 246 ZMAX = MAX( QMAX, ZMAX ) 247* 248* Check for diagonality. 249* 250 IF( E.EQ.ZERO ) THEN 251 DO 20 K = 2, N 252 Z( K ) = Z( 2*K-1 ) 253 20 CONTINUE 254 CALL DLASRT( 'D', N, Z, IINFO ) 255 Z( 2*N-1 ) = D 256 RETURN 257 END IF 258* 259 TRACE = D + E 260* 261* Check for zero data. 262* 263 IF( TRACE.EQ.ZERO ) THEN 264 Z( 2*N-1 ) = ZERO 265 RETURN 266 END IF 267* 268* Check whether the machine is IEEE conformable. 269* 270 IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. 271 $ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 272* 273* Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). 274* 275 DO 30 K = 2*N, 2, -2 276 Z( 2*K ) = ZERO 277 Z( 2*K-1 ) = Z( K ) 278 Z( 2*K-2 ) = ZERO 279 Z( 2*K-3 ) = Z( K-1 ) 280 30 CONTINUE 281* 282 I0 = 1 283 N0 = N 284* 285* Reverse the qd-array, if warranted. 286* 287 IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN 288 IPN4 = 4*( I0+N0 ) 289 DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4 290 TEMP = Z( I4-3 ) 291 Z( I4-3 ) = Z( IPN4-I4-3 ) 292 Z( IPN4-I4-3 ) = TEMP 293 TEMP = Z( I4-1 ) 294 Z( I4-1 ) = Z( IPN4-I4-5 ) 295 Z( IPN4-I4-5 ) = TEMP 296 40 CONTINUE 297 END IF 298* 299* Initial split checking via dqd and Li's test. 300* 301 PP = 0 302* 303 DO 80 K = 1, 2 304* 305 D = Z( 4*N0+PP-3 ) 306 DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 307 IF( Z( I4-1 ).LE.TOL2*D ) THEN 308 Z( I4-1 ) = -ZERO 309 D = Z( I4-3 ) 310 ELSE 311 D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) 312 END IF 313 50 CONTINUE 314* 315* dqd maps Z to ZZ plus Li's test. 316* 317 EMIN = Z( 4*I0+PP+1 ) 318 D = Z( 4*I0+PP-3 ) 319 DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 320 Z( I4-2*PP-2 ) = D + Z( I4-1 ) 321 IF( Z( I4-1 ).LE.TOL2*D ) THEN 322 Z( I4-1 ) = -ZERO 323 Z( I4-2*PP-2 ) = D 324 Z( I4-2*PP ) = ZERO 325 D = Z( I4+1 ) 326 ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. 327 $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN 328 TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) 329 Z( I4-2*PP ) = Z( I4-1 )*TEMP 330 D = D*TEMP 331 ELSE 332 Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) 333 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) 334 END IF 335 EMIN = MIN( EMIN, Z( I4-2*PP ) ) 336 60 CONTINUE 337 Z( 4*N0-PP-2 ) = D 338* 339* Now find qmax. 340* 341 QMAX = Z( 4*I0-PP-2 ) 342 DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 343 QMAX = MAX( QMAX, Z( I4 ) ) 344 70 CONTINUE 345* 346* Prepare for the next iteration on K. 347* 348 PP = 1 - PP 349 80 CONTINUE 350* 351* Initialise variables to pass to DLASQ3. 352* 353 TTYPE = 0 354 DMIN1 = ZERO 355 DMIN2 = ZERO 356 DN = ZERO 357 DN1 = ZERO 358 DN2 = ZERO 359 G = ZERO 360 TAU = ZERO 361* 362 ITER = 2 363 NFAIL = 0 364 NDIV = 2*( N0-I0 ) 365* 366 DO 160 IWHILA = 1, N + 1 367 IF( N0.LT.1 ) 368 $ GO TO 170 369* 370* While array unfinished do 371* 372* E(N0) holds the value of SIGMA when submatrix in I0:N0 373* splits from the rest of the array, but is negated. 374* 375 DESIG = ZERO 376 IF( N0.EQ.N ) THEN 377 SIGMA = ZERO 378 ELSE 379 SIGMA = -Z( 4*N0-1 ) 380 END IF 381 IF( SIGMA.LT.ZERO ) THEN 382 INFO = 1 383 RETURN 384 END IF 385* 386* Find last unreduced submatrix's top index I0, find QMAX and 387* EMIN. Find Gershgorin-type bound if Q's much greater than E's. 388* 389 EMAX = ZERO 390 IF( N0.GT.I0 ) THEN 391 EMIN = ABS( Z( 4*N0-5 ) ) 392 ELSE 393 EMIN = ZERO 394 END IF 395 QMIN = Z( 4*N0-3 ) 396 QMAX = QMIN 397 DO 90 I4 = 4*N0, 8, -4 398 IF( Z( I4-5 ).LE.ZERO ) 399 $ GO TO 100 400 IF( QMIN.GE.FOUR*EMAX ) THEN 401 QMIN = MIN( QMIN, Z( I4-3 ) ) 402 EMAX = MAX( EMAX, Z( I4-5 ) ) 403 END IF 404 QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) 405 EMIN = MIN( EMIN, Z( I4-5 ) ) 406 90 CONTINUE 407 I4 = 4 408* 409 100 CONTINUE 410 I0 = I4 / 4 411 PP = 0 412* 413 IF( N0-I0.GT.1 ) THEN 414 DEE = Z( 4*I0-3 ) 415 DEEMIN = DEE 416 KMIN = I0 417 DO 110 I4 = 4*I0+1, 4*N0-3, 4 418 DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) ) 419 IF( DEE.LE.DEEMIN ) THEN 420 DEEMIN = DEE 421 KMIN = ( I4+3 )/4 422 END IF 423 110 CONTINUE 424 IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 425 $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN 426 IPN4 = 4*( I0+N0 ) 427 PP = 2 428 DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4 429 TEMP = Z( I4-3 ) 430 Z( I4-3 ) = Z( IPN4-I4-3 ) 431 Z( IPN4-I4-3 ) = TEMP 432 TEMP = Z( I4-2 ) 433 Z( I4-2 ) = Z( IPN4-I4-2 ) 434 Z( IPN4-I4-2 ) = TEMP 435 TEMP = Z( I4-1 ) 436 Z( I4-1 ) = Z( IPN4-I4-5 ) 437 Z( IPN4-I4-5 ) = TEMP 438 TEMP = Z( I4 ) 439 Z( I4 ) = Z( IPN4-I4-4 ) 440 Z( IPN4-I4-4 ) = TEMP 441 120 CONTINUE 442 END IF 443 END IF 444* 445* Put -(initial shift) into DMIN. 446* 447 DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) 448* 449* Now I0:N0 is unreduced. 450* PP = 0 for ping, PP = 1 for pong. 451* PP = 2 indicates that flipping was applied to the Z array and 452* and that the tests for deflation upon entry in DLASQ3 453* should not be performed. 454* 455 NBIG = 100*( N0-I0+1 ) 456 DO 140 IWHILB = 1, NBIG 457 IF( I0.GT.N0 ) 458 $ GO TO 150 459* 460* While submatrix unfinished take a good dqds step. 461* 462 CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 463 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 464 $ DN2, G, TAU ) 465* 466 PP = 1 - PP 467* 468* When EMIN is very small check for splits. 469* 470 IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN 471 IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. 472 $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN 473 SPLT = I0 - 1 474 QMAX = Z( 4*I0-3 ) 475 EMIN = Z( 4*I0-1 ) 476 OLDEMN = Z( 4*I0 ) 477 DO 130 I4 = 4*I0, 4*( N0-3 ), 4 478 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. 479 $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN 480 Z( I4-1 ) = -SIGMA 481 SPLT = I4 / 4 482 QMAX = ZERO 483 EMIN = Z( I4+3 ) 484 OLDEMN = Z( I4+4 ) 485 ELSE 486 QMAX = MAX( QMAX, Z( I4+1 ) ) 487 EMIN = MIN( EMIN, Z( I4-1 ) ) 488 OLDEMN = MIN( OLDEMN, Z( I4 ) ) 489 END IF 490 130 CONTINUE 491 Z( 4*N0-1 ) = EMIN 492 Z( 4*N0 ) = OLDEMN 493 I0 = SPLT + 1 494 END IF 495 END IF 496* 497 140 CONTINUE 498* 499 INFO = 2 500* 501* Maximum number of iterations exceeded, restore the shift 502* SIGMA and place the new d's and e's in a qd array. 503* This might need to be done for several blocks 504* 505 I1 = I0 506 N1 = N0 507 145 CONTINUE 508 TEMPQ = Z( 4*I0-3 ) 509 Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA 510 DO K = I0+1, N0 511 TEMPE = Z( 4*K-5 ) 512 Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 )) 513 TEMPQ = Z( 4*K-3 ) 514 Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 ) 515 END DO 516* 517* Prepare to do this on the previous block if there is one 518* 519 IF( I1.GT.1 ) THEN 520 N1 = I1-1 521 DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) ) 522 I1 = I1 - 1 523 END DO 524 SIGMA = -Z(4*N1-1) 525 GO TO 145 526 END IF 527 528 DO K = 1, N 529 Z( 2*K-1 ) = Z( 4*K-3 ) 530* 531* Only the block 1..N0 is unfinished. The rest of the e's 532* must be essentially zero, although sometimes other data 533* has been stored in them. 534* 535 IF( K.LT.N0 ) THEN 536 Z( 2*K ) = Z( 4*K-1 ) 537 ELSE 538 Z( 2*K ) = 0 539 END IF 540 END DO 541 RETURN 542* 543* end IWHILB 544* 545 150 CONTINUE 546* 547 160 CONTINUE 548* 549 INFO = 3 550 RETURN 551* 552* end IWHILA 553* 554 170 CONTINUE 555* 556* Move q's to the front. 557* 558 DO 180 K = 2, N 559 Z( K ) = Z( 4*K-3 ) 560 180 CONTINUE 561* 562* Sort and compute sum of eigenvalues. 563* 564 CALL DLASRT( 'D', N, Z, IINFO ) 565* 566 E = ZERO 567 DO 190 K = N, 1, -1 568 E = E + Z( K ) 569 190 CONTINUE 570* 571* Store trace, sum(eigenvalues) and information on performance. 572* 573 Z( 2*N+1 ) = TRACE 574 Z( 2*N+2 ) = E 575 Z( 2*N+3 ) = DBLE( ITER ) 576 Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) 577 Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) 578 RETURN 579* 580* End of DLASQ2 581* 582 END 583