Return to dlaed9.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DLAED9 used by sstedc. Finds the roots of the secular equation and updates the eigenvectors. Used when the original matrix is dense.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLAED9 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
   22: *                          S, LDS, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
   26: *       DOUBLE PRECISION   RHO
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
   30: *      $                   W( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DLAED9 finds the roots of the secular equation, as defined by the
   40: *> values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
   41: *> appropriate calls to DLAED4 and then stores the new matrix of
   42: *> eigenvectors for use in calculating the next level of Z vectors.
   43: *> \endverbatim
   44: *
   45: *  Arguments:
   46: *  ==========
   47: *
   48: *> \param[in] K
   49: *> \verbatim
   50: *>          K is INTEGER
   51: *>          The number of terms in the rational function to be solved by
   52: *>          DLAED4.  K >= 0.
   53: *> \endverbatim
   54: *>
   55: *> \param[in] KSTART
   56: *> \verbatim
   57: *>          KSTART is INTEGER
   58: *> \endverbatim
   59: *>
   60: *> \param[in] KSTOP
   61: *> \verbatim
   62: *>          KSTOP is INTEGER
   63: *>          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
   64: *>          are to be computed.  1 <= KSTART <= KSTOP <= K.
   65: *> \endverbatim
   66: *>
   67: *> \param[in] N
   68: *> \verbatim
   69: *>          N is INTEGER
   70: *>          The number of rows and columns in the Q matrix.
   71: *>          N >= K (delation may result in N > K).
   72: *> \endverbatim
   73: *>
   74: *> \param[out] D
   75: *> \verbatim
   76: *>          D is DOUBLE PRECISION array, dimension (N)
   77: *>          D(I) contains the updated eigenvalues
   78: *>          for KSTART <= I <= KSTOP.
   79: *> \endverbatim
   80: *>
   81: *> \param[out] Q
   82: *> \verbatim
   83: *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
   84: *> \endverbatim
   85: *>
   86: *> \param[in] LDQ
   87: *> \verbatim
   88: *>          LDQ is INTEGER
   89: *>          The leading dimension of the array Q.  LDQ >= max( 1, N ).
   90: *> \endverbatim
   91: *>
   92: *> \param[in] RHO
   93: *> \verbatim
   94: *>          RHO is DOUBLE PRECISION
   95: *>          The value of the parameter in the rank one update equation.
   96: *>          RHO >= 0 required.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] DLAMDA
  100: *> \verbatim
  101: *>          DLAMDA is DOUBLE PRECISION array, dimension (K)
  102: *>          The first K elements of this array contain the old roots
  103: *>          of the deflated updating problem.  These are the poles
  104: *>          of the secular equation.
  105: *> \endverbatim
  106: *>
  107: *> \param[in] W
  108: *> \verbatim
  109: *>          W is DOUBLE PRECISION array, dimension (K)
  110: *>          The first K elements of this array contain the components
  111: *>          of the deflation-adjusted updating vector.
  112: *> \endverbatim
  113: *>
  114: *> \param[out] S
  115: *> \verbatim
  116: *>          S is DOUBLE PRECISION array, dimension (LDS, K)
  117: *>          Will contain the eigenvectors of the repaired matrix which
  118: *>          will be stored for subsequent Z vector calculation and
  119: *>          multiplied by the previously accumulated eigenvectors
  120: *>          to update the system.
  121: *> \endverbatim
  122: *>
  123: *> \param[in] LDS
  124: *> \verbatim
  125: *>          LDS is INTEGER
  126: *>          The leading dimension of S.  LDS >= max( 1, K ).
  127: *> \endverbatim
  128: *>
  129: *> \param[out] INFO
  130: *> \verbatim
  131: *>          INFO is INTEGER
  132: *>          = 0:  successful exit.
  133: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  134: *>          > 0:  if INFO = 1, an eigenvalue did not converge
  135: *> \endverbatim
  136: *
  137: *  Authors:
  138: *  ========
  139: *
  140: *> \author Univ. of Tennessee 
  141: *> \author Univ. of California Berkeley 
  142: *> \author Univ. of Colorado Denver 
  143: *> \author NAG Ltd. 
  144: *
  145: *> \date September 2012
  146: *
  147: *> \ingroup auxOTHERcomputational
  148: *
  149: *> \par Contributors:
  150: *  ==================
  151: *>
  152: *> Jeff Rutter, Computer Science Division, University of California
  153: *> at Berkeley, USA
  154: *
  155: *  =====================================================================
  156:       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
  157:      $                   S, LDS, INFO )
  158: *
  159: *  -- LAPACK computational routine (version 3.4.2) --
  160: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  161: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  162: *     September 2012
  163: *
  164: *     .. Scalar Arguments ..
  165:       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
  166:       DOUBLE PRECISION   RHO
  167: *     ..
  168: *     .. Array Arguments ..
  169:       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
  170:      $                   W( * )
  171: *     ..
  172: *
  173: *  =====================================================================
  174: *
  175: *     .. Local Scalars ..
  176:       INTEGER            I, J
  177:       DOUBLE PRECISION   TEMP
  178: *     ..
  179: *     .. External Functions ..
  180:       DOUBLE PRECISION   DLAMC3, DNRM2
  181:       EXTERNAL           DLAMC3, DNRM2
  182: *     ..
  183: *     .. External Subroutines ..
  184:       EXTERNAL           DCOPY, DLAED4, XERBLA
  185: *     ..
  186: *     .. Intrinsic Functions ..
  187:       INTRINSIC          MAX, SIGN, SQRT
  188: *     ..
  189: *     .. Executable Statements ..
  190: *
  191: *     Test the input parameters.
  192: *
  193:       INFO = 0
  194: *
  195:       IF( K.LT.0 ) THEN
  196:          INFO = -1
  197:       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
  198:          INFO = -2
  199:       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
  200:      $          THEN
  201:          INFO = -3
  202:       ELSE IF( N.LT.K ) THEN
  203:          INFO = -4
  204:       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
  205:          INFO = -7
  206:       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
  207:          INFO = -12
  208:       END IF
  209:       IF( INFO.NE.0 ) THEN
  210:          CALL XERBLA( 'DLAED9', -INFO )
  211:          RETURN
  212:       END IF
  213: *
  214: *     Quick return if possible
  215: *
  216:       IF( K.EQ.0 )
  217:      $   RETURN
  218: *
  219: *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
  220: *     be computed with high relative accuracy (barring over/underflow).
  221: *     This is a problem on machines without a guard digit in
  222: *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
  223: *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
  224: *     which on any of these machines zeros out the bottommost
  225: *     bit of DLAMDA(I) if it is 1; this makes the subsequent
  226: *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
  227: *     occurs. On binary machines with a guard digit (almost all
  228: *     machines) it does not change DLAMDA(I) at all. On hexadecimal
  229: *     and decimal machines with a guard digit, it slightly
  230: *     changes the bottommost bits of DLAMDA(I). It does not account
  231: *     for hexadecimal or decimal machines without guard digits
  232: *     (we know of none). We use a subroutine call to compute
  233: *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
  234: *     this code.
  235: *
  236:       DO 10 I = 1, N
  237:          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
  238:    10 CONTINUE
  239: *
  240:       DO 20 J = KSTART, KSTOP
  241:          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
  242: *
  243: *        If the zero finder fails, the computation is terminated.
  244: *
  245:          IF( INFO.NE.0 )
  246:      $      GO TO 120
  247:    20 CONTINUE
  248: *
  249:       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
  250:          DO 40 I = 1, K
  251:             DO 30 J = 1, K
  252:                S( J, I ) = Q( J, I )
  253:    30       CONTINUE
  254:    40    CONTINUE
  255:          GO TO 120
  256:       END IF
  257: *
  258: *     Compute updated W.
  259: *
  260:       CALL DCOPY( K, W, 1, S, 1 )
  261: *
  262: *     Initialize W(I) = Q(I,I)
  263: *
  264:       CALL DCOPY( K, Q, LDQ+1, W, 1 )
  265:       DO 70 J = 1, K
  266:          DO 50 I = 1, J - 1
  267:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
  268:    50    CONTINUE
  269:          DO 60 I = J + 1, K
  270:             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
  271:    60    CONTINUE
  272:    70 CONTINUE
  273:       DO 80 I = 1, K
  274:          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
  275:    80 CONTINUE
  276: *
  277: *     Compute eigenvectors of the modified rank-1 modification.
  278: *
  279:       DO 110 J = 1, K
  280:          DO 90 I = 1, K
  281:             Q( I, J ) = W( I ) / Q( I, J )
  282:    90    CONTINUE
  283:          TEMP = DNRM2( K, Q( 1, J ), 1 )
  284:          DO 100 I = 1, K
  285:             S( I, J ) = Q( I, J ) / TEMP
  286:   100    CONTINUE
  287:   110 CONTINUE
  288: *
  289:   120 CONTINUE
  290:       RETURN
  291: *
  292: *     End of DLAED9
  293: *
  294:       END