src/LAPACK/dsterf.f

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