Return to dlaed9.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

    1:       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
    2:      $                   S, LDS, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
   11:       DOUBLE PRECISION   RHO
   12: *     ..
   13: *     .. Array Arguments ..
   14:       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
   15:      $                   W( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DLAED9 finds the roots of the secular equation, as defined by the
   22: *  values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
   23: *  appropriate calls to DLAED4 and then stores the new matrix of
   24: *  eigenvectors for use in calculating the next level of Z vectors.
   25: *
   26: *  Arguments
   27: *  =========
   28: *
   29: *  K       (input) INTEGER
   30: *          The number of terms in the rational function to be solved by
   31: *          DLAED4.  K >= 0.
   32: *
   33: *  KSTART  (input) INTEGER
   34: *  KSTOP   (input) INTEGER
   35: *          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
   36: *          are to be computed.  1 <= KSTART <= KSTOP <= K.
   37: *
   38: *  N       (input) INTEGER
   39: *          The number of rows and columns in the Q matrix.
   40: *          N >= K (delation may result in N > K).
   41: *
   42: *  D       (output) DOUBLE PRECISION array, dimension (N)
   43: *          D(I) contains the updated eigenvalues
   44: *          for KSTART <= I <= KSTOP.
   45: *
   46: *  Q       (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
   47: *
   48: *  LDQ     (input) INTEGER
   49: *          The leading dimension of the array Q.  LDQ >= max( 1, N ).
   50: *
   51: *  RHO     (input) DOUBLE PRECISION
   52: *          The value of the parameter in the rank one update equation.
   53: *          RHO >= 0 required.
   54: *
   55: *  DLAMDA  (input) DOUBLE PRECISION array, dimension (K)
   56: *          The first K elements of this array contain the old roots
   57: *          of the deflated updating problem.  These are the poles
   58: *          of the secular equation.
   59: *
   60: *  W       (input) DOUBLE PRECISION array, dimension (K)
   61: *          The first K elements of this array contain the components
   62: *          of the deflation-adjusted updating vector.
   63: *
   64: *  S       (output) DOUBLE PRECISION array, dimension (LDS, K)
   65: *          Will contain the eigenvectors of the repaired matrix which
   66: *          will be stored for subsequent Z vector calculation and
   67: *          multiplied by the previously accumulated eigenvectors
   68: *          to update the system.
   69: *
   70: *  LDS     (input) INTEGER
   71: *          The leading dimension of S.  LDS >= max( 1, K ).
   72: *
   73: *  INFO    (output) INTEGER
   74: *          = 0:  successful exit.
   75: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
   76: *          > 0:  if INFO = 1, an eigenvalue did not converge
   77: *
   78: *  Further Details
   79: *  ===============
   80: *
   81: *  Based on contributions by
   82: *     Jeff Rutter, Computer Science Division, University of California
   83: *     at Berkeley, USA
   84: *
   85: *  =====================================================================
   86: *
   87: *     .. Local Scalars ..
   88:       INTEGER            I, J
   89:       DOUBLE PRECISION   TEMP
   90: *     ..
   91: *     .. External Functions ..
   92:       DOUBLE PRECISION   DLAMC3, DNRM2
   93:       EXTERNAL           DLAMC3, DNRM2
   94: *     ..
   95: *     .. External Subroutines ..
   96:       EXTERNAL           DCOPY, DLAED4, XERBLA
   97: *     ..
   98: *     .. Intrinsic Functions ..
   99:       INTRINSIC          MAX, SIGN, SQRT
  100: *     ..
  101: *     .. Executable Statements ..
  102: *
  103: *     Test the input parameters.
  104: *
  105:       INFO = 0
  106: *
  107:       IF( K.LT.0 ) THEN
  108:          INFO = -1
  109:       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
  110:          INFO = -2
  111:       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
  112:      $          THEN
  113:          INFO = -3
  114:       ELSE IF( N.LT.K ) THEN
  115:          INFO = -4
  116:       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
  117:          INFO = -7
  118:       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
  119:          INFO = -12
  120:       END IF
  121:       IF( INFO.NE.0 ) THEN
  122:          CALL XERBLA( 'DLAED9', -INFO )
  123:          RETURN
  124:       END IF
  125: *
  126: *     Quick return if possible
  127: *
  128:       IF( K.EQ.0 )
  129:      $   RETURN
  130: *
  131: *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
  132: *     be computed with high relative accuracy (barring over/underflow).
  133: *     This is a problem on machines without a guard digit in
  134: *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
  135: *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
  136: *     which on any of these machines zeros out the bottommost
  137: *     bit of DLAMDA(I) if it is 1; this makes the subsequent
  138: *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
  139: *     occurs. On binary machines with a guard digit (almost all
  140: *     machines) it does not change DLAMDA(I) at all. On hexadecimal
  141: *     and decimal machines with a guard digit, it slightly
  142: *     changes the bottommost bits of DLAMDA(I). It does not account
  143: *     for hexadecimal or decimal machines without guard digits
  144: *     (we know of none). We use a subroutine call to compute
  145: *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
  146: *     this code.
  147: *
  148:       DO 10 I = 1, N
  149:          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
  150:    10 CONTINUE
  151: *
  152:       DO 20 J = KSTART, KSTOP
  153:          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
  154: *
  155: *        If the zero finder fails, the computation is terminated.
  156: *
  157:          IF( INFO.NE.0 )
  158:      $      GO TO 120
  159:    20 CONTINUE
  160: *
  161:       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
  162:          DO 40 I = 1, K
  163:             DO 30 J = 1, K
  164:                S( J, I ) = Q( J, I )
  165:    30       CONTINUE
  166:    40    CONTINUE
  167:          GO TO 120
  168:       END IF
  169: *
  170: *     Compute updated W.
  171: *
  172:       CALL DCOPY( K, W, 1, S, 1 )
  173: *
  174: *     Initialize W(I) = Q(I,I)
  175: *
  176:       CALL DCOPY( K, Q, LDQ+1, W, 1 )
  177:       DO 70 J = 1, K
  178:          DO 50 I = 1, J - 1
  179:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
  180:    50    CONTINUE
  181:          DO 60 I = J + 1, K
  182:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
  183:    60    CONTINUE
  184:    70 CONTINUE
  185:       DO 80 I = 1, K
  186:          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
  187:    80 CONTINUE
  188: *
  189: *     Compute eigenvectors of the modified rank-1 modification.
  190: *
  191:       DO 110 J = 1, K
  192:          DO 90 I = 1, K
  193:             Q( I, J ) = W( I ) / Q( I, J )
  194:    90    CONTINUE
  195:          TEMP = DNRM2( K, Q( 1, J ), 1 )
  196:          DO 100 I = 1, K
  197:             S( I, J ) = Q( I, J ) / TEMP
  198:   100    CONTINUE
  199:   110 CONTINUE
  200: *
  201:   120 CONTINUE
  202:       RETURN
  203: *
  204: *     End of DLAED9
  205: *
  206:       END