Return to zlaed7.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b ZLAED7 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. 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 ZLAED7 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed7.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed7.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed7.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
   22: *                          LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
   23: *                          GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
   24: *                          INFO )
   25: * 
   26: *       .. Scalar Arguments ..
   27: *       INTEGER            CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
   28: *      $                   TLVLS
   29: *       DOUBLE PRECISION   RHO
   30: *       ..
   31: *       .. Array Arguments ..
   32: *       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
   33: *      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
   34: *       DOUBLE PRECISION   D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
   35: *       COMPLEX*16         Q( LDQ, * ), WORK( * )
   36: *       ..
   37: *  
   38: *
   39: *> \par Purpose:
   40: *  =============
   41: *>
   42: *> \verbatim
   43: *>
   44: *> ZLAED7 computes the updated eigensystem of a diagonal
   45: *> matrix after modification by a rank-one symmetric matrix. This
   46: *> routine is used only for the eigenproblem which requires all
   47: *> eigenvalues and optionally eigenvectors of a dense or banded
   48: *> Hermitian matrix that has been reduced to tridiagonal form.
   49: *>
   50: *>   T = Q(in) ( D(in) + RHO * Z*Z**H ) Q**H(in) = Q(out) * D(out) * Q**H(out)
   51: *>
   52: *>   where Z = Q**Hu, u is a vector of length N with ones in the
   53: *>   CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
   54: *>
   55: *>    The eigenvectors of the original matrix are stored in Q, and the
   56: *>    eigenvalues are in D.  The algorithm consists of three stages:
   57: *>
   58: *>       The first stage consists of deflating the size of the problem
   59: *>       when there are multiple eigenvalues or if there is a zero in
   60: *>       the Z vector.  For each such occurence the dimension of the
   61: *>       secular equation problem is reduced by one.  This stage is
   62: *>       performed by the routine DLAED2.
   63: *>
   64: *>       The second stage consists of calculating the updated
   65: *>       eigenvalues. This is done by finding the roots of the secular
   66: *>       equation via the routine DLAED4 (as called by SLAED3).
   67: *>       This routine also calculates the eigenvectors of the current
   68: *>       problem.
   69: *>
   70: *>       The final stage consists of computing the updated eigenvectors
   71: *>       directly using the updated eigenvalues.  The eigenvectors for
   72: *>       the current problem are multiplied with the eigenvectors from
   73: *>       the overall problem.
   74: *> \endverbatim
   75: *
   76: *  Arguments:
   77: *  ==========
   78: *
   79: *> \param[in] N
   80: *> \verbatim
   81: *>          N is INTEGER
   82: *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
   83: *> \endverbatim
   84: *>
   85: *> \param[in] CUTPNT
   86: *> \verbatim
   87: *>          CUTPNT is INTEGER
   88: *>         Contains the location of the last eigenvalue in the leading
   89: *>         sub-matrix.  min(1,N) <= CUTPNT <= N.
   90: *> \endverbatim
   91: *>
   92: *> \param[in] QSIZ
   93: *> \verbatim
   94: *>          QSIZ is INTEGER
   95: *>         The dimension of the unitary matrix used to reduce
   96: *>         the full matrix to tridiagonal form.  QSIZ >= N.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] TLVLS
  100: *> \verbatim
  101: *>          TLVLS is INTEGER
  102: *>         The total number of merging levels in the overall divide and
  103: *>         conquer tree.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] CURLVL
  107: *> \verbatim
  108: *>          CURLVL is INTEGER
  109: *>         The current level in the overall merge routine,
  110: *>         0 <= curlvl <= tlvls.
  111: *> \endverbatim
  112: *>
  113: *> \param[in] CURPBM
  114: *> \verbatim
  115: *>          CURPBM is INTEGER
  116: *>         The current problem in the current level in the overall
  117: *>         merge routine (counting from upper left to lower right).
  118: *> \endverbatim
  119: *>
  120: *> \param[in,out] D
  121: *> \verbatim
  122: *>          D is DOUBLE PRECISION array, dimension (N)
  123: *>         On entry, the eigenvalues of the rank-1-perturbed matrix.
  124: *>         On exit, the eigenvalues of the repaired matrix.
  125: *> \endverbatim
  126: *>
  127: *> \param[in,out] Q
  128: *> \verbatim
  129: *>          Q is COMPLEX*16 array, dimension (LDQ,N)
  130: *>         On entry, the eigenvectors of the rank-1-perturbed matrix.
  131: *>         On exit, the eigenvectors of the repaired tridiagonal matrix.
  132: *> \endverbatim
  133: *>
  134: *> \param[in] LDQ
  135: *> \verbatim
  136: *>          LDQ is INTEGER
  137: *>         The leading dimension of the array Q.  LDQ >= max(1,N).
  138: *> \endverbatim
  139: *>
  140: *> \param[in] RHO
  141: *> \verbatim
  142: *>          RHO is DOUBLE PRECISION
  143: *>         Contains the subdiagonal element used to create the rank-1
  144: *>         modification.
  145: *> \endverbatim
  146: *>
  147: *> \param[out] INDXQ
  148: *> \verbatim
  149: *>          INDXQ is INTEGER array, dimension (N)
  150: *>         This contains the permutation which will reintegrate the
  151: *>         subproblem just solved back into sorted order,
  152: *>         ie. D( INDXQ( I = 1, N ) ) will be in ascending order.
  153: *> \endverbatim
  154: *>
  155: *> \param[out] IWORK
  156: *> \verbatim
  157: *>          IWORK is INTEGER array, dimension (4*N)
  158: *> \endverbatim
  159: *>
  160: *> \param[out] RWORK
  161: *> \verbatim
  162: *>          RWORK is DOUBLE PRECISION array,
  163: *>                                 dimension (3*N+2*QSIZ*N)
  164: *> \endverbatim
  165: *>
  166: *> \param[out] WORK
  167: *> \verbatim
  168: *>          WORK is COMPLEX*16 array, dimension (QSIZ*N)
  169: *> \endverbatim
  170: *>
  171: *> \param[in,out] QSTORE
  172: *> \verbatim
  173: *>          QSTORE is DOUBLE PRECISION array, dimension (N**2+1)
  174: *>         Stores eigenvectors of submatrices encountered during
  175: *>         divide and conquer, packed together. QPTR points to
  176: *>         beginning of the submatrices.
  177: *> \endverbatim
  178: *>
  179: *> \param[in,out] QPTR
  180: *> \verbatim
  181: *>          QPTR is INTEGER array, dimension (N+2)
  182: *>         List of indices pointing to beginning of submatrices stored
  183: *>         in QSTORE. The submatrices are numbered starting at the
  184: *>         bottom left of the divide and conquer tree, from left to
  185: *>         right and bottom to top.
  186: *> \endverbatim
  187: *>
  188: *> \param[in] PRMPTR
  189: *> \verbatim
  190: *>          PRMPTR is INTEGER array, dimension (N lg N)
  191: *>         Contains a list of pointers which indicate where in PERM a
  192: *>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
  193: *>         indicates the size of the permutation and also the size of
  194: *>         the full, non-deflated problem.
  195: *> \endverbatim
  196: *>
  197: *> \param[in] PERM
  198: *> \verbatim
  199: *>          PERM is INTEGER array, dimension (N lg N)
  200: *>         Contains the permutations (from deflation and sorting) to be
  201: *>         applied to each eigenblock.
  202: *> \endverbatim
  203: *>
  204: *> \param[in] GIVPTR
  205: *> \verbatim
  206: *>          GIVPTR is INTEGER array, dimension (N lg N)
  207: *>         Contains a list of pointers which indicate where in GIVCOL a
  208: *>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
  209: *>         indicates the number of Givens rotations.
  210: *> \endverbatim
  211: *>
  212: *> \param[in] GIVCOL
  213: *> \verbatim
  214: *>          GIVCOL is INTEGER array, dimension (2, N lg N)
  215: *>         Each pair of numbers indicates a pair of columns to take place
  216: *>         in a Givens rotation.
  217: *> \endverbatim
  218: *>
  219: *> \param[in] GIVNUM
  220: *> \verbatim
  221: *>          GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
  222: *>         Each number indicates the S value to be used in the
  223: *>         corresponding Givens rotation.
  224: *> \endverbatim
  225: *>
  226: *> \param[out] INFO
  227: *> \verbatim
  228: *>          INFO is INTEGER
  229: *>          = 0:  successful exit.
  230: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  231: *>          > 0:  if INFO = 1, an eigenvalue did not converge
  232: *> \endverbatim
  233: *
  234: *  Authors:
  235: *  ========
  236: *
  237: *> \author Univ. of Tennessee 
  238: *> \author Univ. of California Berkeley 
  239: *> \author Univ. of Colorado Denver 
  240: *> \author NAG Ltd. 
  241: *
  242: *> \date September 2012
  243: *
  244: *> \ingroup complex16OTHERcomputational
  245: *
  246: *  =====================================================================
  247:       SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
  248:      $                   LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
  249:      $                   GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
  250:      $                   INFO )
  251: *
  252: *  -- LAPACK computational routine (version 3.4.2) --
  253: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  254: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  255: *     September 2012
  256: *
  257: *     .. Scalar Arguments ..
  258:       INTEGER            CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
  259:      $                   TLVLS
  260:       DOUBLE PRECISION   RHO
  261: *     ..
  262: *     .. Array Arguments ..
  263:       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
  264:      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
  265:       DOUBLE PRECISION   D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
  266:       COMPLEX*16         Q( LDQ, * ), WORK( * )
  267: *     ..
  268: *
  269: *  =====================================================================
  270: *
  271: *     .. Local Scalars ..
  272:       INTEGER            COLTYP, CURR, I, IDLMDA, INDX,
  273:      $                   INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
  274: *     ..
  275: *     .. External Subroutines ..
  276:       EXTERNAL           DLAED9, DLAEDA, DLAMRG, XERBLA, ZLACRM, ZLAED8
  277: *     ..
  278: *     .. Intrinsic Functions ..
  279:       INTRINSIC          MAX, MIN
  280: *     ..
  281: *     .. Executable Statements ..
  282: *
  283: *     Test the input parameters.
  284: *
  285:       INFO = 0
  286: *
  287: *     IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
  288: *        INFO = -1
  289: *     ELSE IF( N.LT.0 ) THEN
  290:       IF( N.LT.0 ) THEN
  291:          INFO = -1
  292:       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
  293:          INFO = -2
  294:       ELSE IF( QSIZ.LT.N ) THEN
  295:          INFO = -3
  296:       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
  297:          INFO = -9
  298:       END IF
  299:       IF( INFO.NE.0 ) THEN
  300:          CALL XERBLA( 'ZLAED7', -INFO )
  301:          RETURN
  302:       END IF
  303: *
  304: *     Quick return if possible
  305: *
  306:       IF( N.EQ.0 )
  307:      $   RETURN
  308: *
  309: *     The following values are for bookkeeping purposes only.  They are
  310: *     integer pointers which indicate the portion of the workspace
  311: *     used by a particular array in DLAED2 and SLAED3.
  312: *
  313:       IZ = 1
  314:       IDLMDA = IZ + N
  315:       IW = IDLMDA + N
  316:       IQ = IW + N
  317: *
  318:       INDX = 1
  319:       INDXC = INDX + N
  320:       COLTYP = INDXC + N
  321:       INDXP = COLTYP + N
  322: *
  323: *     Form the z-vector which consists of the last row of Q_1 and the
  324: *     first row of Q_2.
  325: *
  326:       PTR = 1 + 2**TLVLS
  327:       DO 10 I = 1, CURLVL - 1
  328:          PTR = PTR + 2**( TLVLS-I )
  329:    10 CONTINUE
  330:       CURR = PTR + CURPBM
  331:       CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
  332:      $             GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
  333:      $             RWORK( IZ+N ), INFO )
  334: *
  335: *     When solving the final problem, we no longer need the stored data,
  336: *     so we will overwrite the data from this level onto the previously
  337: *     used storage space.
  338: *
  339:       IF( CURLVL.EQ.TLVLS ) THEN
  340:          QPTR( CURR ) = 1
  341:          PRMPTR( CURR ) = 1
  342:          GIVPTR( CURR ) = 1
  343:       END IF
  344: *
  345: *     Sort and Deflate eigenvalues.
  346: *
  347:       CALL ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
  348:      $             RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
  349:      $             IWORK( INDXP ), IWORK( INDX ), INDXQ,
  350:      $             PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
  351:      $             GIVCOL( 1, GIVPTR( CURR ) ),
  352:      $             GIVNUM( 1, GIVPTR( CURR ) ), INFO )
  353:       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
  354:       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
  355: *
  356: *     Solve Secular Equation.
  357: *
  358:       IF( K.NE.0 ) THEN
  359:          CALL DLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
  360:      $                RWORK( IDLMDA ), RWORK( IW ),
  361:      $                QSTORE( QPTR( CURR ) ), K, INFO )
  362:          CALL ZLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
  363:      $                LDQ, RWORK( IQ ) )
  364:          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
  365:          IF( INFO.NE.0 ) THEN
  366:             RETURN
  367:          END IF
  368: *
  369: *     Prepare the INDXQ sorting premutation.
  370: *
  371:          N1 = K
  372:          N2 = N - K
  373:          CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
  374:       ELSE
  375:          QPTR( CURR+1 ) = QPTR( CURR )
  376:          DO 20 I = 1, N
  377:             INDXQ( I ) = I
  378:    20    CONTINUE
  379:       END IF
  380: *
  381:       RETURN
  382: *
  383: *     End of ZLAED7
  384: *
  385:       END