vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlasqtest/dlasq2.f

   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