Return to dsterf.f CVS log

Up to [local] / rpl / lapack / lapack

Commit initial.

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