github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dlasqtest/dlasq2.f (about) 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 September 2012 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.4.2) -- 116 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 117 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 118 * September 2012 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 291 TEMP = Z( I4-3 ) 292 Z( I4-3 ) = Z( IPN4-I4-3 ) 293 Z( IPN4-I4-3 ) = TEMP 294 TEMP = Z( I4-1 ) 295 Z( I4-1 ) = Z( IPN4-I4-5 ) 296 Z( IPN4-I4-5 ) = TEMP 297 40 CONTINUE 298 END IF 299 * 300 * Initial split checking via dqd and Li's test. 301 * 302 PP = 0 303 * 304 DO 80 K = 1, 2 305 * 306 D = Z( 4*N0+PP-3 ) 307 DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 308 IF( Z( I4-1 ).LE.TOL2*D ) THEN 309 Z( I4-1 ) = -ZERO 310 D = Z( I4-3 ) 311 ELSE 312 D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) 313 END IF 314 50 CONTINUE 315 * 316 * dqd maps Z to ZZ plus Li's test. 317 * 318 EMIN = Z( 4*I0+PP+1 ) 319 D = Z( 4*I0+PP-3 ) 320 DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 321 Z( I4-2*PP-2 ) = D + Z( I4-1 ) 322 IF( Z( I4-1 ).LE.TOL2*D ) THEN 323 Z( I4-1 ) = -ZERO 324 Z( I4-2*PP-2 ) = D 325 Z( I4-2*PP ) = ZERO 326 D = Z( I4+1 ) 327 ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. 328 $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN 329 TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) 330 Z( I4-2*PP ) = Z( I4-1 )*TEMP 331 D = D*TEMP 332 ELSE 333 Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) 334 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) 335 END IF 336 EMIN = MIN( EMIN, Z( I4-2*PP ) ) 337 60 CONTINUE 338 Z( 4*N0-PP-2 ) = D 339 * 340 * Now find qmax. 341 * 342 QMAX = Z( 4*I0-PP-2 ) 343 DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 344 QMAX = MAX( QMAX, Z( I4 ) ) 345 70 CONTINUE 346 * 347 * Prepare for the next iteration on K. 348 * 349 PP = 1 - PP 350 80 CONTINUE 351 352 * 353 * Initialise variables to pass to DLASQ3. 354 * 355 TTYPE = 0 356 DMIN1 = ZERO 357 DMIN2 = ZERO 358 DN = ZERO 359 DN1 = ZERO 360 DN2 = ZERO 361 G = ZERO 362 TAU = ZERO 363 * 364 ITER = 2 365 NFAIL = 0 366 NDIV = 2*( N0-I0 ) 367 * 368 DO 160 IWHILA = 1, N + 1 369 370 IF( N0.LT.1 ) THEN 371 GO TO 170 372 END IF 373 * 374 * While array unfinished do 375 * 376 * E(N0) holds the value of SIGMA when submatrix in I0:N0 377 * splits from the rest of the array, but is negated. 378 * 379 DESIG = ZERO 380 IF( N0.EQ.N ) THEN 381 SIGMA = ZERO 382 ELSE 383 SIGMA = -Z( 4*N0-1 ) 384 END IF 385 IF( SIGMA.LT.ZERO ) THEN 386 INFO = 1 387 RETURN 388 END IF 389 * 390 * Find last unreduced submatrix's top index I0, find QMAX and 391 * EMIN. Find Gershgorin-type bound if Q's much greater than E's. 392 * 393 EMAX = ZERO 394 IF( N0.GT.I0 ) THEN 395 EMIN = ABS( Z( 4*N0-5 ) ) 396 ELSE 397 EMIN = ZERO 398 END IF 399 QMIN = Z( 4*N0-3 ) 400 QMAX = QMIN 401 DO 90 I4 = 4*N0, 8, -4 402 IF( Z( I4-5 ).LE.ZERO ) 403 $ GO TO 100 404 IF( QMIN.GE.FOUR*EMAX ) THEN 405 QMIN = MIN( QMIN, Z( I4-3 ) ) 406 EMAX = MAX( EMAX, Z( I4-5 ) ) 407 END IF 408 QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) 409 EMIN = MIN( EMIN, Z( I4-5 ) ) 410 90 CONTINUE 411 I4 = 4 412 * 413 100 CONTINUE 414 I0 = I4 / 4 415 416 PP = 0 417 * 418 IF( N0-I0.GT.1 ) THEN 419 DEE = Z( 4*I0-3 ) 420 DEEMIN = DEE 421 KMIN = I0 422 DO 110 I4 = 4*I0+1, 4*N0-3, 4 423 DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) ) 424 IF( DEE.LE.DEEMIN ) THEN 425 DEEMIN = DEE 426 KMIN = ( I4+3 )/4 427 END IF 428 110 CONTINUE 429 IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 430 $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN 431 IPN4 = 4*( I0+N0 ) 432 PP = 2 433 DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4 434 TEMP = Z( I4-3 ) 435 Z( I4-3 ) = Z( IPN4-I4-3 ) 436 Z( IPN4-I4-3 ) = TEMP 437 TEMP = Z( I4-2 ) 438 Z( I4-2 ) = Z( IPN4-I4-2 ) 439 Z( IPN4-I4-2 ) = TEMP 440 TEMP = Z( I4-1 ) 441 Z( I4-1 ) = Z( IPN4-I4-5 ) 442 Z( IPN4-I4-5 ) = TEMP 443 TEMP = Z( I4 ) 444 Z( I4 ) = Z( IPN4-I4-4 ) 445 Z( IPN4-I4-4 ) = TEMP 446 120 CONTINUE 447 END IF 448 END IF 449 * 450 * Put -(initial shift) into DMIN. 451 * 452 DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) 453 * 454 * Now I0:N0 is unreduced. 455 * PP = 0 for ping, PP = 1 for pong. 456 * PP = 2 indicates that flipping was applied to the Z array and 457 * and that the tests for deflation upon entry in DLASQ3 458 * should not be performed. 459 * 460 NBIG = 100*( N0-I0+1 ) 461 DO 140 IWHILB = 1, NBIG 462 463 IF( I0.GT.N0 ) 464 $ GO TO 150 465 * 466 467 ! Print out test cases 468 469 write(3,*) "{" 470 write(3,*) "i0: ", I0, "," 471 write(3,*) "n0: ", N0, "," 472 write(3,'(9999(g0))',advance="no") "z: []float64{" 473 do i = 1, 4*n 474 write (3,'(99999(e24.16,a))',advance="no") z(i), "," 475 end do 476 write (3,*) "}," 477 write (3,*) "pp: ", PP, "," 478 write (3,*) "dmin: ", DMIN, "," 479 write (3,*) "desig:", DESIG, "," 480 write (3,*) "qmax: ", QMAX, "," 481 write (3,*) "ttype:", TTYPE, "," 482 write (3,*) "dmin1:", DMIN1, "," 483 write (3,*) "dmin2:", DMIN2, "," 484 write (3,*) "dn: ", DN, "," 485 write (3,*) "dn1: ", DN1, "," 486 write (3,*) "dn2: ", DN2, "," 487 write (3,*) "g: ", G, "," 488 write (3,*) "tau: ", TAU, "," 489 write (3,*) "nFail:", NFAIL, "," 490 write (3,*) "iter: ", ITER, "," 491 write (3,*) "sigma:", SIGMA, "," 492 write (3,*) "nDiv: ", NDIV, "," 493 494 * While submatrix unfinished take a good dqds step. 495 * 496 497 498 CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 499 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 500 $ DN2, G, TAU ) 501 502 503 ! Write the outputs 504 write(3,'(9999(g0))',advance="no") "zOut: []float64{" 505 do i = 1, 4*n 506 write (3,'(99999(e24.16,a))',advance="no") z(i), "," 507 end do 508 write (3,*) "}," 509 write (3,*) "i0Out:",I0, "," 510 write (3,*) "n0Out:", N0, "," 511 write (3,*) "ppOut:", PP, "," 512 write (3,*) "dminOut:", DMIN, "," 513 write (3,*) "desigOut:", DESIG, "," 514 write (3,*) "sigmaOut:", SIGMA, "," 515 write (3,*) "qmaxOut:", QMAX, "," 516 write (3,*) "nFailOut:", NFAIL, "," 517 write (3,*) "iterOut:", ITER, "," 518 write (3,*) "nDivOut:", NDIV, "," 519 write (3,*) "ttypeOut:", TTYPE, "," 520 write (3,*) "dmin1Out:", DMIN1, "," 521 write (3,*) "dmin2Out:", DMIN2, "," 522 write (3,*) "dnOut:", DN, "," 523 write (3,*) "dn1Out:", DN1, "," 524 write (3,*) "dn2Out:", DN2, "," 525 write (3,*) "gOut:", G, "," 526 write (3,*) "tauOut:", TAU, "," 527 528 write (3,*) "}," 529 530 531 PP = 1 - PP 532 * 533 * When EMIN is very small check for splits. 534 * 535 IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN 536 IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. 537 $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN 538 SPLT = I0 - 1 539 QMAX = Z( 4*I0-3 ) 540 EMIN = Z( 4*I0-1 ) 541 OLDEMN = Z( 4*I0 ) 542 DO 130 I4 = 4*I0, 4*( N0-3 ), 4 543 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. 544 $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN 545 Z( I4-1 ) = -SIGMA 546 SPLT = I4 / 4 547 QMAX = ZERO 548 EMIN = Z( I4+3 ) 549 OLDEMN = Z( I4+4 ) 550 ELSE 551 QMAX = MAX( QMAX, Z( I4+1 ) ) 552 EMIN = MIN( EMIN, Z( I4-1 ) ) 553 OLDEMN = MIN( OLDEMN, Z( I4 ) ) 554 END IF 555 130 CONTINUE 556 Z( 4*N0-1 ) = EMIN 557 Z( 4*N0 ) = OLDEMN 558 I0 = SPLT + 1 559 END IF 560 END IF 561 * 562 140 CONTINUE 563 * 564 INFO = 2 565 * 566 * Maximum number of iterations exceeded, restore the shift 567 * SIGMA and place the new d's and e's in a qd array. 568 * This might need to be done for several blocks 569 * 570 I1 = I0 571 N1 = N0 572 145 CONTINUE 573 574 TEMPQ = Z( 4*I0-3 ) 575 Z( 4*I0-3 ) = Z( 4*I0-3 ) + SIGMA 576 DO K = I0+1, N0 577 TEMPE = Z( 4*K-5 ) 578 Z( 4*K-5 ) = Z( 4*K-5 ) * (TEMPQ / Z( 4*K-7 )) 579 TEMPQ = Z( 4*K-3 ) 580 Z( 4*K-3 ) = Z( 4*K-3 ) + SIGMA + TEMPE - Z( 4*K-5 ) 581 END DO 582 * 583 * Prepare to do this on the previous block if there is one 584 * 585 IF( I1.GT.1 ) THEN 586 N1 = I1-1 587 DO WHILE( ( I1.GE.2 ) .AND. ( Z(4*I1-5).GE.ZERO ) ) 588 I1 = I1 - 1 589 END DO 590 SIGMA = -Z(4*N1-1) 591 GO TO 145 592 END IF 593 594 DO K = 1, N 595 Z( 2*K-1 ) = Z( 4*K-3 ) 596 * 597 * Only the block 1..N0 is unfinished. The rest of the e's 598 * must be essentially zero, although sometimes other data 599 * has been stored in them. 600 * 601 IF( K.LT.N0 ) THEN 602 Z( 2*K ) = Z( 4*K-1 ) 603 ELSE 604 Z( 2*K ) = 0 605 END IF 606 END DO 607 RETURN 608 * 609 * end IWHILB 610 * 611 150 CONTINUE 612 * 613 160 CONTINUE 614 * 615 INFO = 3 616 RETURN 617 * 618 * end IWHILA 619 * 620 170 CONTINUE 621 * 622 623 * Move q's to the front. 624 * 625 DO 180 K = 2, N 626 Z( K ) = Z( 4*K-3 ) 627 180 CONTINUE 628 * 629 * Sort and compute sum of eigenvalues. 630 * 631 CALL DLASRT( 'D', N, Z, IINFO ) 632 * 633 634 E = ZERO 635 DO 190 K = N, 1, -1 636 E = E + Z( K ) 637 190 CONTINUE 638 * 639 * Store trace, sum(eigenvalues) and information on performance. 640 * 641 642 Z( 2*N+1 ) = TRACE 643 Z( 2*N+2 ) = E 644 Z( 2*N+3 ) = DBLE( ITER ) 645 Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) 646 Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) 647 648 RETURN 649 * 650 * End of DLASQ2 651 * 652 END