Return to zlaed7.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
    2:      $                   LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
    3:      $                   GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
    4:      $                   INFO )
    5: *
    6: *  -- LAPACK routine (version 3.2) --
    7: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    8: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    9: *     November 2006
   10: *
   11: *     .. Scalar Arguments ..
   12:       INTEGER            CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
   13:      $                   TLVLS
   14:       DOUBLE PRECISION   RHO
   15: *     ..
   16: *     .. Array Arguments ..
   17:       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
   18:      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
   19:       DOUBLE PRECISION   D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
   20:       COMPLEX*16         Q( LDQ, * ), WORK( * )
   21: *     ..
   22: *
   23: *  Purpose
   24: *  =======
   25: *
   26: *  ZLAED7 computes the updated eigensystem of a diagonal
   27: *  matrix after modification by a rank-one symmetric matrix. This
   28: *  routine is used only for the eigenproblem which requires all
   29: *  eigenvalues and optionally eigenvectors of a dense or banded
   30: *  Hermitian matrix that has been reduced to tridiagonal form.
   31: *
   32: *    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
   33: *
   34: *    where Z = Q'u, u is a vector of length N with ones in the
   35: *    CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
   36: *
   37: *     The eigenvectors of the original matrix are stored in Q, and the
   38: *     eigenvalues are in D.  The algorithm consists of three stages:
   39: *
   40: *        The first stage consists of deflating the size of the problem
   41: *        when there are multiple eigenvalues or if there is a zero in
   42: *        the Z vector.  For each such occurence the dimension of the
   43: *        secular equation problem is reduced by one.  This stage is
   44: *        performed by the routine DLAED2.
   45: *
   46: *        The second stage consists of calculating the updated
   47: *        eigenvalues. This is done by finding the roots of the secular
   48: *        equation via the routine DLAED4 (as called by SLAED3).
   49: *        This routine also calculates the eigenvectors of the current
   50: *        problem.
   51: *
   52: *        The final stage consists of computing the updated eigenvectors
   53: *        directly using the updated eigenvalues.  The eigenvectors for
   54: *        the current problem are multiplied with the eigenvectors from
   55: *        the overall problem.
   56: *
   57: *  Arguments
   58: *  =========
   59: *
   60: *  N      (input) INTEGER
   61: *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
   62: *
   63: *  CUTPNT (input) INTEGER
   64: *         Contains the location of the last eigenvalue in the leading
   65: *         sub-matrix.  min(1,N) <= CUTPNT <= N.
   66: *
   67: *  QSIZ   (input) INTEGER
   68: *         The dimension of the unitary matrix used to reduce
   69: *         the full matrix to tridiagonal form.  QSIZ >= N.
   70: *
   71: *  TLVLS  (input) INTEGER
   72: *         The total number of merging levels in the overall divide and
   73: *         conquer tree.
   74: *
   75: *  CURLVL (input) INTEGER
   76: *         The current level in the overall merge routine,
   77: *         0 <= curlvl <= tlvls.
   78: *
   79: *  CURPBM (input) INTEGER
   80: *         The current problem in the current level in the overall
   81: *         merge routine (counting from upper left to lower right).
   82: *
   83: *  D      (input/output) DOUBLE PRECISION array, dimension (N)
   84: *         On entry, the eigenvalues of the rank-1-perturbed matrix.
   85: *         On exit, the eigenvalues of the repaired matrix.
   86: *
   87: *  Q      (input/output) COMPLEX*16 array, dimension (LDQ,N)
   88: *         On entry, the eigenvectors of the rank-1-perturbed matrix.
   89: *         On exit, the eigenvectors of the repaired tridiagonal matrix.
   90: *
   91: *  LDQ    (input) INTEGER
   92: *         The leading dimension of the array Q.  LDQ >= max(1,N).
   93: *
   94: *  RHO    (input) DOUBLE PRECISION
   95: *         Contains the subdiagonal element used to create the rank-1
   96: *         modification.
   97: *
   98: *  INDXQ  (output) INTEGER array, dimension (N)
   99: *         This contains the permutation which will reintegrate the
  100: *         subproblem just solved back into sorted order,
  101: *         ie. D( INDXQ( I = 1, N ) ) will be in ascending order.
  102: *
  103: *  IWORK  (workspace) INTEGER array, dimension (4*N)
  104: *
  105: *  RWORK  (workspace) DOUBLE PRECISION array,
  106: *                                 dimension (3*N+2*QSIZ*N)
  107: *
  108: *  WORK   (workspace) COMPLEX*16 array, dimension (QSIZ*N)
  109: *
  110: *  QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1)
  111: *         Stores eigenvectors of submatrices encountered during
  112: *         divide and conquer, packed together. QPTR points to
  113: *         beginning of the submatrices.
  114: *
  115: *  QPTR   (input/output) INTEGER array, dimension (N+2)
  116: *         List of indices pointing to beginning of submatrices stored
  117: *         in QSTORE. The submatrices are numbered starting at the
  118: *         bottom left of the divide and conquer tree, from left to
  119: *         right and bottom to top.
  120: *
  121: *  PRMPTR (input) INTEGER array, dimension (N lg N)
  122: *         Contains a list of pointers which indicate where in PERM a
  123: *         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
  124: *         indicates the size of the permutation and also the size of
  125: *         the full, non-deflated problem.
  126: *
  127: *  PERM   (input) INTEGER array, dimension (N lg N)
  128: *         Contains the permutations (from deflation and sorting) to be
  129: *         applied to each eigenblock.
  130: *
  131: *  GIVPTR (input) INTEGER array, dimension (N lg N)
  132: *         Contains a list of pointers which indicate where in GIVCOL a
  133: *         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
  134: *         indicates the number of Givens rotations.
  135: *
  136: *  GIVCOL (input) INTEGER array, dimension (2, N lg N)
  137: *         Each pair of numbers indicates a pair of columns to take place
  138: *         in a Givens rotation.
  139: *
  140: *  GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
  141: *         Each number indicates the S value to be used in the
  142: *         corresponding Givens rotation.
  143: *
  144: *  INFO   (output) INTEGER
  145: *          = 0:  successful exit.
  146: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  147: *          > 0:  if INFO = 1, an eigenvalue did not converge
  148: *
  149: *  =====================================================================
  150: *
  151: *     .. Local Scalars ..
  152:       INTEGER            COLTYP, CURR, I, IDLMDA, INDX,
  153:      $                   INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
  154: *     ..
  155: *     .. External Subroutines ..
  156:       EXTERNAL           DLAED9, DLAEDA, DLAMRG, XERBLA, ZLACRM, ZLAED8
  157: *     ..
  158: *     .. Intrinsic Functions ..
  159:       INTRINSIC          MAX, MIN
  160: *     ..
  161: *     .. Executable Statements ..
  162: *
  163: *     Test the input parameters.
  164: *
  165:       INFO = 0
  166: *
  167: *     IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
  168: *        INFO = -1
  169: *     ELSE IF( N.LT.0 ) THEN
  170:       IF( N.LT.0 ) THEN
  171:          INFO = -1
  172:       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
  173:          INFO = -2
  174:       ELSE IF( QSIZ.LT.N ) THEN
  175:          INFO = -3
  176:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
  177:          INFO = -9
  178:       END IF
  179:       IF( INFO.NE.0 ) THEN
  180:          CALL XERBLA( 'ZLAED7', -INFO )
  181:          RETURN
  182:       END IF
  183: *
  184: *     Quick return if possible
  185: *
  186:       IF( N.EQ.0 )
  187:      $   RETURN
  188: *
  189: *     The following values are for bookkeeping purposes only.  They are
  190: *     integer pointers which indicate the portion of the workspace
  191: *     used by a particular array in DLAED2 and SLAED3.
  192: *
  193:       IZ = 1
  194:       IDLMDA = IZ + N
  195:       IW = IDLMDA + N
  196:       IQ = IW + N
  197: *
  198:       INDX = 1
  199:       INDXC = INDX + N
  200:       COLTYP = INDXC + N
  201:       INDXP = COLTYP + N
  202: *
  203: *     Form the z-vector which consists of the last row of Q_1 and the
  204: *     first row of Q_2.
  205: *
  206:       PTR = 1 + 2**TLVLS
  207:       DO 10 I = 1, CURLVL - 1
  208:          PTR = PTR + 2**( TLVLS-I )
  209:    10 CONTINUE
  210:       CURR = PTR + CURPBM
  211:       CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
  212:      $             GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
  213:      $             RWORK( IZ+N ), INFO )
  214: *
  215: *     When solving the final problem, we no longer need the stored data,
  216: *     so we will overwrite the data from this level onto the previously
  217: *     used storage space.
  218: *
  219:       IF( CURLVL.EQ.TLVLS ) THEN
  220:          QPTR( CURR ) = 1
  221:          PRMPTR( CURR ) = 1
  222:          GIVPTR( CURR ) = 1
  223:       END IF
  224: *
  225: *     Sort and Deflate eigenvalues.
  226: *
  227:       CALL ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
  228:      $             RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
  229:      $             IWORK( INDXP ), IWORK( INDX ), INDXQ,
  230:      $             PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
  231:      $             GIVCOL( 1, GIVPTR( CURR ) ),
  232:      $             GIVNUM( 1, GIVPTR( CURR ) ), INFO )
  233:       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
  234:       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
  235: *
  236: *     Solve Secular Equation.
  237: *
  238:       IF( K.NE.0 ) THEN
  239:          CALL DLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
  240:      $                RWORK( IDLMDA ), RWORK( IW ),
  241:      $                QSTORE( QPTR( CURR ) ), K, INFO )
  242:          CALL ZLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
  243:      $                LDQ, RWORK( IQ ) )
  244:          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
  245:          IF( INFO.NE.0 ) THEN
  246:             RETURN
  247:          END IF
  248: *
  249: *     Prepare the INDXQ sorting premutation.
  250: *
  251:          N1 = K
  252:          N2 = N - K
  253:          CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
  254:       ELSE
  255:          QPTR( CURR+1 ) = QPTR( CURR )
  256:          DO 20 I = 1, N
  257:             INDXQ( I ) = I
  258:    20    CONTINUE
  259:       END IF
  260: *
  261:       RETURN
  262: *
  263: *     End of ZLAED7
  264: *
  265:       END