SRC/dsterf.f

   1 *> \brief \b DSTERF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSTERF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSTERF( N, D, E, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION   D( * ), E( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DSTERF 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTERF( 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       DOUBLE PRECISION   D( * ), E( * )
  99 *     ..
 100 *
 101 *  =====================================================================
 102 *
 103 *     .. Parameters ..
 104       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
 105       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 106      $                   THREE = 3.0D0 )
 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       DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
 114      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
 115      $                   SIGMA, SSFMAX, SSFMIN, RMAX
 116 *     ..
 117 *     .. External Functions ..
 118       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
 119       EXTERNAL           DLAMCH, DLANST, DLAPY2
 120 *     ..
 121 *     .. External Subroutines ..
 122       EXTERNAL           DLAE2, DLASCL, DLASRT, 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( 'DSTERF', -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 = DLAMCH( 'E' )
 146       EPS2 = EPS**2
 147       SAFMIN = DLAMCH( 'S' )
 148       SAFMAX = ONE / SAFMIN
 149       SSFMAX = SQRT( SAFMAX ) / THREE
 150       SSFMIN = SQRT( SAFMIN ) / EPS2
 151       RMAX = DLAMCH( 'O' )
 152 *
 153 *     Compute the eigenvalues of the tridiagonal matrix.
 154 *
 155       NMAXIT = N*MAXIT
 156       SIGMA = ZERO
 157       JTOT = 0
 158 *
 159 *     Determine where the matrix splits and choose QL or QR iteration
 160 *     for each block, according to whether top or bottom diagonal
 161 *     element is smaller.
 162 *
 163       L1 = 1
 164 *
 165    10 CONTINUE
 166       IF( L1.GT.N )
 167      $   GO TO 170
 168       IF( L1.GT.1 )
 169      $   E( L1-1 ) = ZERO
 170       DO 20 M = L1, N - 1
 171          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
 172      $       1 ) ) ) )*EPS ) THEN
 173             E( M ) = ZERO
 174             GO TO 30
 175          END IF
 176    20 CONTINUE
 177       M = N
 178 *
 179    30 CONTINUE
 180       L = L1
 181       LSV = L
 182       LEND = M
 183       LENDSV = LEND
 184       L1 = M + 1
 185       IF( LEND.EQ.L )
 186      $   GO TO 10
 187 *
 188 *     Scale submatrix in rows and columns L to LEND
 189 *
 190       ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
 191       ISCALE = 0
 192       IF( ANORM.EQ.ZERO )
 193      $   GO TO 10
 194       IF( (ANORM.GT.SSFMAX) ) THEN
 195          ISCALE = 1
 196          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
 197      $                INFO )
 198          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
 199      $                INFO )
 200       ELSE IF( ANORM.LT.SSFMIN ) THEN
 201          ISCALE = 2
 202          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
 203      $                INFO )
 204          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
 205      $                INFO )
 206       END IF
 207 *
 208       DO 40 I = L, LEND - 1
 209          E( I ) = E( I )**2
 210    40 CONTINUE
 211 *
 212 *     Choose between QL and QR iteration
 213 *
 214       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
 215          LEND = LSV
 216          L = LENDSV
 217       END IF
 218 *
 219       IF( LEND.GE.L ) THEN
 220 *
 221 *        QL Iteration
 222 *
 223 *        Look for small subdiagonal element.
 224 *
 225    50    CONTINUE
 226          IF( L.NE.LEND ) THEN
 227             DO 60 M = L, LEND - 1
 228                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
 229      $            GO TO 70
 230    60       CONTINUE
 231          END IF
 232          M = LEND
 233 *
 234    70    CONTINUE
 235          IF( M.LT.LEND )
 236      $      E( M ) = ZERO
 237          P = D( L )
 238          IF( M.EQ.L )
 239      $      GO TO 90
 240 *
 241 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
 242 *        eigenvalues.
 243 *
 244          IF( M.EQ.L+1 ) THEN
 245             RTE = SQRT( E( L ) )
 246             CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
 247             D( L ) = RT1
 248             D( L+1 ) = RT2
 249             E( L ) = ZERO
 250             L = L + 2
 251             IF( L.LE.LEND )
 252      $         GO TO 50
 253             GO TO 150
 254          END IF
 255 *
 256          IF( JTOT.EQ.NMAXIT )
 257      $      GO TO 150
 258          JTOT = JTOT + 1
 259 *
 260 *        Form shift.
 261 *
 262          RTE = SQRT( E( L ) )
 263          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
 264          R = DLAPY2( SIGMA, ONE )
 265          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 266 *
 267          C = ONE
 268          S = ZERO
 269          GAMMA = D( M ) - SIGMA
 270          P = GAMMA*GAMMA
 271 *
 272 *        Inner loop
 273 *
 274          DO 80 I = M - 1, L, -1
 275             BB = E( I )
 276             R = P + BB
 277             IF( I.NE.M-1 )
 278      $         E( I+1 ) = S*R
 279             OLDC = C
 280             C = P / R
 281             S = BB / R
 282             OLDGAM = GAMMA
 283             ALPHA = D( I )
 284             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 285             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
 286             IF( C.NE.ZERO ) THEN
 287                P = ( GAMMA*GAMMA ) / C
 288             ELSE
 289                P = OLDC*BB
 290             END IF
 291    80    CONTINUE
 292 *
 293          E( L ) = S*P
 294          D( L ) = SIGMA + GAMMA
 295          GO TO 50
 296 *
 297 *        Eigenvalue found.
 298 *
 299    90    CONTINUE
 300          D( L ) = P
 301 *
 302          L = L + 1
 303          IF( L.LE.LEND )
 304      $      GO TO 50
 305          GO TO 150
 306 *
 307       ELSE
 308 *
 309 *        QR Iteration
 310 *
 311 *        Look for small superdiagonal element.
 312 *
 313   100    CONTINUE
 314          DO 110 M = L, LEND + 1, -1
 315             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
 316      $         GO TO 120
 317   110    CONTINUE
 318          M = LEND
 319 *
 320   120    CONTINUE
 321          IF( M.GT.LEND )
 322      $      E( M-1 ) = ZERO
 323          P = D( L )
 324          IF( M.EQ.L )
 325      $      GO TO 140
 326 *
 327 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
 328 *        eigenvalues.
 329 *
 330          IF( M.EQ.L-1 ) THEN
 331             RTE = SQRT( E( L-1 ) )
 332             CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
 333             D( L ) = RT1
 334             D( L-1 ) = RT2
 335             E( L-1 ) = ZERO
 336             L = L - 2
 337             IF( L.GE.LEND )
 338      $         GO TO 100
 339             GO TO 150
 340          END IF
 341 *
 342          IF( JTOT.EQ.NMAXIT )
 343      $      GO TO 150
 344          JTOT = JTOT + 1
 345 *
 346 *        Form shift.
 347 *
 348          RTE = SQRT( E( L-1 ) )
 349          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
 350          R = DLAPY2( SIGMA, ONE )
 351          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 352 *
 353          C = ONE
 354          S = ZERO
 355          GAMMA = D( M ) - SIGMA
 356          P = GAMMA*GAMMA
 357 *
 358 *        Inner loop
 359 *
 360          DO 130 I = M, L - 1
 361             BB = E( I )
 362             R = P + BB
 363             IF( I.NE.M )
 364      $         E( I-1 ) = S*R
 365             OLDC = C
 366             C = P / R
 367             S = BB / R
 368             OLDGAM = GAMMA
 369             ALPHA = D( I+1 )
 370             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 371             D( I ) = OLDGAM + ( ALPHA-GAMMA )
 372             IF( C.NE.ZERO ) THEN
 373                P = ( GAMMA*GAMMA ) / C
 374             ELSE
 375                P = OLDC*BB
 376             END IF
 377   130    CONTINUE
 378 *
 379          E( L-1 ) = S*P
 380          D( L ) = SIGMA + GAMMA
 381          GO TO 100
 382 *
 383 *        Eigenvalue found.
 384 *
 385   140    CONTINUE
 386          D( L ) = P
 387 *
 388          L = L - 1
 389          IF( L.GE.LEND )
 390      $      GO TO 100
 391          GO TO 150
 392 *
 393       END IF
 394 *
 395 *     Undo scaling if necessary
 396 *
 397   150 CONTINUE
 398       IF( ISCALE.EQ.1 )
 399      $   CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
 400      $                D( LSV ), N, INFO )
 401       IF( ISCALE.EQ.2 )
 402      $   CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
 403      $                D( LSV ), N, INFO )
 404 *
 405 *     Check for no convergence to an eigenvalue after a total
 406 *     of N*MAXIT iterations.
 407 *
 408       IF( JTOT.LT.NMAXIT )
 409      $   GO TO 10
 410       DO 160 I = 1, N - 1
 411          IF( E( I ).NE.ZERO )
 412      $      INFO = INFO + 1
 413   160 CONTINUE
 414       GO TO 180
 415 *
 416 *     Sort eigenvalues in increasing order.
 417 *
 418   170 CONTINUE
 419       CALL DLASRT( 'I', N, D, INFO )
 420 *
 421   180 CONTINUE
 422       RETURN
 423 *
 424 *     End of DSTERF
 425 *
 426       END