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