LAPACK 3.3.0

claed7.f

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