SRC/ssterf.f

   1 *> \brief \b SSTERF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSTERF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssterf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssterf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssterf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSTERF( N, D, E, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       REAL               D( * ), E( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> SSTERF computes all eigenvalues of a symmetric tridiagonal matrix
  37 *> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] N
  44 *> \verbatim
  45 *>          N is INTEGER
  46 *>          The order of the matrix.  N >= 0.
  47 *> \endverbatim
  48 *>
  49 *> \param[in,out] D
  50 *> \verbatim
  51 *>          D is REAL array, dimension (N)
  52 *>          On entry, the n diagonal elements of the tridiagonal matrix.
  53 *>          On exit, if INFO = 0, the eigenvalues in ascending order.
  54 *> \endverbatim
  55 *>
  56 *> \param[in,out] E
  57 *> \verbatim
  58 *>          E is REAL array, dimension (N-1)
  59 *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
  60 *>          matrix.
  61 *>          On exit, E has been destroyed.
  62 *> \endverbatim
  63 *>
  64 *> \param[out] INFO
  65 *> \verbatim
  66 *>          INFO is INTEGER
  67 *>          = 0:  successful exit
  68 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  69 *>          > 0:  the algorithm failed to find all of the eigenvalues in
  70 *>                a total of 30*N iterations; if INFO = i, then i
  71 *>                elements of E have not converged to zero.
  72 *> \endverbatim
  73 *
  74 *  Authors:
  75 *  ========
  76 *
  77 *> \author Univ. of Tennessee
  78 *> \author Univ. of California Berkeley
  79 *> \author Univ. of Colorado Denver
  80 *> \author NAG Ltd.
  81 *
  82 *> \date November 2011
  83 *
  84 *> \ingroup auxOTHERcomputational
  85 *
  86 *  =====================================================================
  87       SUBROUTINE SSTERF( N, D, E, INFO )
  88 *
  89 *  -- LAPACK computational routine (version 3.4.0) --
  90 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  91 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  92 *     November 2011
  93 *
  94 *     .. Scalar Arguments ..
  95       INTEGER            INFO, N
  96 *     ..
  97 *     .. Array Arguments ..
  98       REAL               D( * ), E( * )
  99 *     ..
 100 *
 101 *  =====================================================================
 102 *
 103 *     .. Parameters ..
 104       REAL               ZERO, ONE, TWO, THREE
 105       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
 106      $                   THREE = 3.0E0 )
 107       INTEGER            MAXIT
 108       PARAMETER          ( MAXIT = 30 )
 109 *     ..
 110 *     .. Local Scalars ..
 111       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
 112      $                   NMAXIT
 113       REAL               ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
 114      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
 115      $                   SIGMA, SSFMAX, SSFMIN
 116 *     ..
 117 *     .. External Functions ..
 118       REAL               SLAMCH, SLANST, SLAPY2
 119       EXTERNAL           SLAMCH, SLANST, SLAPY2
 120 *     ..
 121 *     .. External Subroutines ..
 122       EXTERNAL           SLAE2, SLASCL, SLASRT, XERBLA
 123 *     ..
 124 *     .. Intrinsic Functions ..
 125       INTRINSIC          ABS, SIGN, SQRT
 126 *     ..
 127 *     .. Executable Statements ..
 128 *
 129 *     Test the input parameters.
 130 *
 131       INFO = 0
 132 *
 133 *     Quick return if possible
 134 *
 135       IF( N.LT.0 ) THEN
 136          INFO = -1
 137          CALL XERBLA( 'SSTERF', -INFO )
 138          RETURN
 139       END IF
 140       IF( N.LE.1 )
 141      $   RETURN
 142 *
 143 *     Determine the unit roundoff for this environment.
 144 *
 145       EPS = SLAMCH( 'E' )
 146       EPS2 = EPS**2
 147       SAFMIN = SLAMCH( 'S' )
 148       SAFMAX = ONE / SAFMIN
 149       SSFMAX = SQRT( SAFMAX ) / THREE
 150       SSFMIN = SQRT( SAFMIN ) / EPS2
 151 *
 152 *     Compute the eigenvalues of the tridiagonal matrix.
 153 *
 154       NMAXIT = N*MAXIT
 155       SIGMA = ZERO
 156       JTOT = 0
 157 *
 158 *     Determine where the matrix splits and choose QL or QR iteration
 159 *     for each block, according to whether top or bottom diagonal
 160 *     element is smaller.
 161 *
 162       L1 = 1
 163 *
 164    10 CONTINUE
 165       IF( L1.GT.N )
 166      $   GO TO 170
 167       IF( L1.GT.1 )
 168      $   E( L1-1 ) = ZERO
 169       DO 20 M = L1, N - 1
 170          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*
 171      $       SQRT( ABS( D( M+1 ) ) ) )*EPS ) THEN
 172             E( M ) = ZERO
 173             GO TO 30
 174          END IF
 175    20 CONTINUE
 176       M = N
 177 *
 178    30 CONTINUE
 179       L = L1
 180       LSV = L
 181       LEND = M
 182       LENDSV = LEND
 183       L1 = M + 1
 184       IF( LEND.EQ.L )
 185      $   GO TO 10
 186 *
 187 *     Scale submatrix in rows and columns L to LEND
 188 *
 189       ANORM = SLANST( 'M', LEND-L+1, D( L ), E( L ) )
 190       ISCALE = 0
 191       IF( ANORM.EQ.ZERO )
 192      $   GO TO 10
 193       IF( ANORM.GT.SSFMAX ) THEN
 194          ISCALE = 1
 195          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
 196      $                INFO )
 197          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
 198      $                INFO )
 199       ELSE IF( ANORM.LT.SSFMIN ) THEN
 200          ISCALE = 2
 201          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
 202      $                INFO )
 203          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
 204      $                INFO )
 205       END IF
 206 *
 207       DO 40 I = L, LEND - 1
 208          E( I ) = E( I )**2
 209    40 CONTINUE
 210 *
 211 *     Choose between QL and QR iteration
 212 *
 213       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
 214          LEND = LSV
 215          L = LENDSV
 216       END IF
 217 *
 218       IF( LEND.GE.L ) THEN
 219 *
 220 *        QL Iteration
 221 *
 222 *        Look for small subdiagonal element.
 223 *
 224    50    CONTINUE
 225          IF( L.NE.LEND ) THEN
 226             DO 60 M = L, LEND - 1
 227                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
 228      $            GO TO 70
 229    60       CONTINUE
 230          END IF
 231          M = LEND
 232 *
 233    70    CONTINUE
 234          IF( M.LT.LEND )
 235      $      E( M ) = ZERO
 236          P = D( L )
 237          IF( M.EQ.L )
 238      $      GO TO 90
 239 *
 240 *        If remaining matrix is 2 by 2, use SLAE2 to compute its
 241 *        eigenvalues.
 242 *
 243          IF( M.EQ.L+1 ) THEN
 244             RTE = SQRT( E( L ) )
 245             CALL SLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
 246             D( L ) = RT1
 247             D( L+1 ) = RT2
 248             E( L ) = ZERO
 249             L = L + 2
 250             IF( L.LE.LEND )
 251      $         GO TO 50
 252             GO TO 150
 253          END IF
 254 *
 255          IF( JTOT.EQ.NMAXIT )
 256      $      GO TO 150
 257          JTOT = JTOT + 1
 258 *
 259 *        Form shift.
 260 *
 261          RTE = SQRT( E( L ) )
 262          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
 263          R = SLAPY2( SIGMA, ONE )
 264          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 265 *
 266          C = ONE
 267          S = ZERO
 268          GAMMA = D( M ) - SIGMA
 269          P = GAMMA*GAMMA
 270 *
 271 *        Inner loop
 272 *
 273          DO 80 I = M - 1, L, -1
 274             BB = E( I )
 275             R = P + BB
 276             IF( I.NE.M-1 )
 277      $         E( I+1 ) = S*R
 278             OLDC = C
 279             C = P / R
 280             S = BB / R
 281             OLDGAM = GAMMA
 282             ALPHA = D( I )
 283             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 284             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
 285             IF( C.NE.ZERO ) THEN
 286                P = ( GAMMA*GAMMA ) / C
 287             ELSE
 288                P = OLDC*BB
 289             END IF
 290    80    CONTINUE
 291 *
 292          E( L ) = S*P
 293          D( L ) = SIGMA + GAMMA
 294          GO TO 50
 295 *
 296 *        Eigenvalue found.
 297 *
 298    90    CONTINUE
 299          D( L ) = P
 300 *
 301          L = L + 1
 302          IF( L.LE.LEND )
 303      $      GO TO 50
 304          GO TO 150
 305 *
 306       ELSE
 307 *
 308 *        QR Iteration
 309 *
 310 *        Look for small superdiagonal element.
 311 *
 312   100    CONTINUE
 313          DO 110 M = L, LEND + 1, -1
 314             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
 315      $         GO TO 120
 316   110    CONTINUE
 317          M = LEND
 318 *
 319   120    CONTINUE
 320          IF( M.GT.LEND )
 321      $      E( M-1 ) = ZERO
 322          P = D( L )
 323          IF( M.EQ.L )
 324      $      GO TO 140
 325 *
 326 *        If remaining matrix is 2 by 2, use SLAE2 to compute its
 327 *        eigenvalues.
 328 *
 329          IF( M.EQ.L-1 ) THEN
 330             RTE = SQRT( E( L-1 ) )
 331             CALL SLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
 332             D( L ) = RT1
 333             D( L-1 ) = RT2
 334             E( L-1 ) = ZERO
 335             L = L - 2
 336             IF( L.GE.LEND )
 337      $         GO TO 100
 338             GO TO 150
 339          END IF
 340 *
 341          IF( JTOT.EQ.NMAXIT )
 342      $      GO TO 150
 343          JTOT = JTOT + 1
 344 *
 345 *        Form shift.
 346 *
 347          RTE = SQRT( E( L-1 ) )
 348          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
 349          R = SLAPY2( SIGMA, ONE )
 350          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 351 *
 352          C = ONE
 353          S = ZERO
 354          GAMMA = D( M ) - SIGMA
 355          P = GAMMA*GAMMA
 356 *
 357 *        Inner loop
 358 *
 359          DO 130 I = M, L - 1
 360             BB = E( I )
 361             R = P + BB
 362             IF( I.NE.M )
 363      $         E( I-1 ) = S*R
 364             OLDC = C
 365             C = P / R
 366             S = BB / R
 367             OLDGAM = GAMMA
 368             ALPHA = D( I+1 )
 369             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 370             D( I ) = OLDGAM + ( ALPHA-GAMMA )
 371             IF( C.NE.ZERO ) THEN
 372                P = ( GAMMA*GAMMA ) / C
 373             ELSE
 374                P = OLDC*BB
 375             END IF
 376   130    CONTINUE
 377 *
 378          E( L-1 ) = S*P
 379          D( L ) = SIGMA + GAMMA
 380          GO TO 100
 381 *
 382 *        Eigenvalue found.
 383 *
 384   140    CONTINUE
 385          D( L ) = P
 386 *
 387          L = L - 1
 388          IF( L.GE.LEND )
 389      $      GO TO 100
 390          GO TO 150
 391 *
 392       END IF
 393 *
 394 *     Undo scaling if necessary
 395 *
 396   150 CONTINUE
 397       IF( ISCALE.EQ.1 )
 398      $   CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
 399      $                D( LSV ), N, INFO )
 400       IF( ISCALE.EQ.2 )
 401      $   CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
 402      $                D( LSV ), N, INFO )
 403 *
 404 *     Check for no convergence to an eigenvalue after a total
 405 *     of N*MAXIT iterations.
 406 *
 407       IF( JTOT.LT.NMAXIT )
 408      $   GO TO 10
 409       DO 160 I = 1, N - 1
 410          IF( E( I ).NE.ZERO )
 411      $      INFO = INFO + 1
 412   160 CONTINUE
 413       GO TO 180
 414 *
 415 *     Sort eigenvalues in increasing order.
 416 *
 417   170 CONTINUE
 418       CALL SLASRT( 'I', N, D, INFO )
 419 *
 420   180 CONTINUE
 421       RETURN
 422 *
 423 *     End of SSTERF
 424 *
 425       END