d0/de7/dlaed7_8f_source.html

00001       SUBROUTINE DLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
00002      $                   LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
00003      $                   PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
00004      $                   INFO )
00005 *
00006 *  -- LAPACK routine (version 3.2) --
00007 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00008 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00009 *     November 2006
00010 *
00011 *     .. Scalar Arguments ..
00012       INTEGER            CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
00013      $                   QSIZ, TLVLS
00014       DOUBLE PRECISION   RHO
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
00018      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
00019       DOUBLE PRECISION   D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
00020      $                   QSTORE( * ), WORK( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  DLAED7 computes the updated eigensystem of a diagonal
00027 *  matrix after modification by a rank-one symmetric matrix. This
00028 *  routine is used only for the eigenproblem which requires all
00029 *  eigenvalues and optionally eigenvectors of a dense symmetric matrix
00030 *  that has been reduced to tridiagonal form.  DLAED1 handles
00031 *  the case in which all eigenvalues and eigenvectors of a symmetric
00032 *  tridiagonal matrix are desired.
00033 *
00034 *    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
00035 *
00036 *     where Z = Q'u, u is a vector of length N with ones in the
00037 *     CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
00038 *
00039 *     The eigenvectors of the original matrix are stored in Q, and the
00040 *     eigenvalues are in D.  The algorithm consists of three stages:
00041 *
00042 *        The first stage consists of deflating the size of the problem
00043 *        when there are multiple eigenvalues or if there is a zero in
00044 *        the Z vector.  For each such occurence the dimension of the
00045 *        secular equation problem is reduced by one.  This stage is
00046 *        performed by the routine DLAED8.
00047 *
00048 *        The second stage consists of calculating the updated
00049 *        eigenvalues. This is done by finding the roots of the secular
00050 *        equation via the routine DLAED4 (as called by DLAED9).
00051 *        This routine also calculates the eigenvectors of the current
00052 *        problem.
00053 *
00054 *        The final stage consists of computing the updated eigenvectors
00055 *        directly using the updated eigenvalues.  The eigenvectors for
00056 *        the current problem are multiplied with the eigenvectors from
00057 *        the overall problem.
00058 *
00059 *  Arguments
00060 *  =========
00061 *
00062 *  ICOMPQ  (input) INTEGER
00063 *          = 0:  Compute eigenvalues only.
00064 *          = 1:  Compute eigenvectors of original dense symmetric matrix
00065 *                also.  On entry, Q contains the orthogonal matrix used
00066 *                to reduce the original matrix to tridiagonal form.
00067 *
00068 *  N      (input) INTEGER
00069 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00070 *
00071 *  QSIZ   (input) INTEGER
00072 *         The dimension of the orthogonal matrix used to reduce
00073 *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00074 *
00075 *  TLVLS  (input) INTEGER
00076 *         The total number of merging levels in the overall divide and
00077 *         conquer tree.
00078 *
00079 *  CURLVL (input) INTEGER
00080 *         The current level in the overall merge routine,
00081 *         0 <= CURLVL <= TLVLS.
00082 *
00083 *  CURPBM (input) INTEGER
00084 *         The current problem in the current level in the overall
00085 *         merge routine (counting from upper left to lower right).
00086 *
00087 *  D      (input/output) DOUBLE PRECISION array, dimension (N)
00088 *         On entry, the eigenvalues of the rank-1-perturbed matrix.
00089 *         On exit, the eigenvalues of the repaired matrix.
00090 *
00091 *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
00092 *         On entry, the eigenvectors of the rank-1-perturbed matrix.
00093 *         On exit, the eigenvectors of the repaired tridiagonal matrix.
00094 *
00095 *  LDQ    (input) INTEGER
00096 *         The leading dimension of the array Q.  LDQ >= max(1,N).
00097 *
00098 *  INDXQ  (output) INTEGER array, dimension (N)
00099 *         The permutation which will reintegrate the subproblem just
00100 *         solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
00101 *         will be in ascending order.
00102 *
00103 *  RHO    (input) DOUBLE PRECISION
00104 *         The subdiagonal element used to create the rank-1
00105 *         modification.
00106 *
00107 *  CUTPNT (input) INTEGER
00108 *         Contains the location of the last eigenvalue in the leading
00109 *         sub-matrix.  min(1,N) <= CUTPNT <= N.
00110 *
00111 *  QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1)
00112 *         Stores eigenvectors of submatrices encountered during
00113 *         divide and conquer, packed together. QPTR points to
00114 *         beginning of the submatrices.
00115 *
00116 *  QPTR   (input/output) INTEGER array, dimension (N+2)
00117 *         List of indices pointing to beginning of submatrices stored
00118 *         in QSTORE. The submatrices are numbered starting at the
00119 *         bottom left of the divide and conquer tree, from left to
00120 *         right and bottom to top.
00121 *
00122 *  PRMPTR (input) INTEGER array, dimension (N lg N)
00123 *         Contains a list of pointers which indicate where in PERM a
00124 *         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
00125 *         indicates the size of the permutation and also the size of
00126 *         the full, non-deflated problem.
00127 *
00128 *  PERM   (input) INTEGER array, dimension (N lg N)
00129 *         Contains the permutations (from deflation and sorting) to be
00130 *         applied to each eigenblock.
00131 *
00132 *  GIVPTR (input) INTEGER array, dimension (N lg N)
00133 *         Contains a list of pointers which indicate where in GIVCOL a
00134 *         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
00135 *         indicates the number of Givens rotations.
00136 *
00137 *  GIVCOL (input) INTEGER array, dimension (2, N lg N)
00138 *         Each pair of numbers indicates a pair of columns to take place
00139 *         in a Givens rotation.
00140 *
00141 *  GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
00142 *         Each number indicates the S value to be used in the
00143 *         corresponding Givens rotation.
00144 *
00145 *  WORK   (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N)
00146 *
00147 *  IWORK  (workspace) INTEGER array, dimension (4*N)
00148 *
00149 *  INFO   (output) INTEGER
00150 *          = 0:  successful exit.
00151 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00152 *          > 0:  if INFO = 1, an eigenvalue did not converge
00153 *
00154 *  Further Details
00155 *  ===============
00156 *
00157 *  Based on contributions by
00158 *     Jeff Rutter, Computer Science Division, University of California
00159 *     at Berkeley, USA
00160 *
00161 *  =====================================================================
00162 *
00163 *     .. Parameters ..
00164       DOUBLE PRECISION   ONE, ZERO
00165       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
00166 *     ..
00167 *     .. Local Scalars ..
00168       INTEGER            COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP,
00169      $                   IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR
00170 *     ..
00171 *     .. External Subroutines ..
00172       EXTERNAL           DGEMM, DLAED8, DLAED9, DLAEDA, DLAMRG, XERBLA
00173 *     ..
00174 *     .. Intrinsic Functions ..
00175       INTRINSIC          MAX, MIN
00176 *     ..
00177 *     .. Executable Statements ..
00178 *
00179 *     Test the input parameters.
00180 *
00181       INFO = 0
00182 *
00183       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
00184          INFO = -1
00185       ELSE IF( N.LT.0 ) THEN
00186          INFO = -2
00187       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
00188          INFO = -4
00189       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00190          INFO = -9
00191       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
00192          INFO = -12
00193       END IF
00194       IF( INFO.NE.0 ) THEN
00195          CALL XERBLA( 'DLAED7', -INFO )
00196          RETURN
00197       END IF
00198 *
00199 *     Quick return if possible
00200 *
00201       IF( N.EQ.0 )
00202      $   RETURN
00203 *
00204 *     The following values are for bookkeeping purposes only.  They are
00205 *     integer pointers which indicate the portion of the workspace
00206 *     used by a particular array in DLAED8 and DLAED9.
00207 *
00208       IF( ICOMPQ.EQ.1 ) THEN
00209          LDQ2 = QSIZ
00210       ELSE
00211          LDQ2 = N
00212       END IF
00213 *
00214       IZ = 1
00215       IDLMDA = IZ + N
00216       IW = IDLMDA + N
00217       IQ2 = IW + N
00218       IS = IQ2 + N*LDQ2
00219 *
00220       INDX = 1
00221       INDXC = INDX + N
00222       COLTYP = INDXC + N
00223       INDXP = COLTYP + N
00224 *
00225 *     Form the z-vector which consists of the last row of Q_1 and the
00226 *     first row of Q_2.
00227 *
00228       PTR = 1 + 2**TLVLS
00229       DO 10 I = 1, CURLVL - 1
00230          PTR = PTR + 2**( TLVLS-I )
00231    10 CONTINUE
00232       CURR = PTR + CURPBM
00233       CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
00234      $             GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ),
00235      $             WORK( IZ+N ), INFO )
00236 *
00237 *     When solving the final problem, we no longer need the stored data,
00238 *     so we will overwrite the data from this level onto the previously
00239 *     used storage space.
00240 *
00241       IF( CURLVL.EQ.TLVLS ) THEN
00242          QPTR( CURR ) = 1
00243          PRMPTR( CURR ) = 1
00244          GIVPTR( CURR ) = 1
00245       END IF
00246 *
00247 *     Sort and Deflate eigenvalues.
00248 *
00249       CALL DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT,
00250      $             WORK( IZ ), WORK( IDLMDA ), WORK( IQ2 ), LDQ2,
00251      $             WORK( IW ), PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
00252      $             GIVCOL( 1, GIVPTR( CURR ) ),
00253      $             GIVNUM( 1, GIVPTR( CURR ) ), IWORK( INDXP ),
00254      $             IWORK( INDX ), INFO )
00255       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
00256       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
00257 *
00258 *     Solve Secular Equation.
00259 *
00260       IF( K.NE.0 ) THEN
00261          CALL DLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ),
00262      $                WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO )
00263          IF( INFO.NE.0 )
00264      $      GO TO 30
00265          IF( ICOMPQ.EQ.1 ) THEN
00266             CALL DGEMM( 'N', 'N', QSIZ, K, K, ONE, WORK( IQ2 ), LDQ2,
00267      $                  QSTORE( QPTR( CURR ) ), K, ZERO, Q, LDQ )
00268          END IF
00269          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
00270 *
00271 *     Prepare the INDXQ sorting permutation.
00272 *
00273          N1 = K
00274          N2 = N - K
00275          CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
00276       ELSE
00277          QPTR( CURR+1 ) = QPTR( CURR )
00278          DO 20 I = 1, N
00279             INDXQ( I ) = I
00280    20    CONTINUE
00281       END IF
00282 *
00283    30 CONTINUE
00284       RETURN
00285 *
00286 *     End of DLAED7
00287 *
00288       END