LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slaed7.f
Go to the documentation of this file.
00001 *> \brief \b SLAED7
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLAED7 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed7.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed7.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed7.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
00022 *                          LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
00023 *                          PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
00024 *                          INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       INTEGER            CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
00028 *      $                   QSIZ, TLVLS
00029 *       REAL               RHO
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
00033 *      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
00034 *       REAL               D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
00035 *      $                   QSTORE( * ), WORK( * )
00036 *       ..
00037 *  
00038 *
00039 *> \par Purpose:
00040 *  =============
00041 *>
00042 *> \verbatim
00043 *>
00044 *> SLAED7 computes the updated eigensystem of a diagonal
00045 *> matrix after modification by a rank-one symmetric matrix. This
00046 *> routine is used only for the eigenproblem which requires all
00047 *> eigenvalues and optionally eigenvectors of a dense symmetric matrix
00048 *> that has been reduced to tridiagonal form.  SLAED1 handles
00049 *> the case in which all eigenvalues and eigenvectors of a symmetric
00050 *> tridiagonal matrix are desired.
00051 *>
00052 *>   T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out)
00053 *>
00054 *>    where Z = Q**Tu, u is a vector of length N with ones in the
00055 *>    CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
00056 *>
00057 *>    The eigenvectors of the original matrix are stored in Q, and the
00058 *>    eigenvalues are in D.  The algorithm consists of three stages:
00059 *>
00060 *>       The first stage consists of deflating the size of the problem
00061 *>       when there are multiple eigenvalues or if there is a zero in
00062 *>       the Z vector.  For each such occurence the dimension of the
00063 *>       secular equation problem is reduced by one.  This stage is
00064 *>       performed by the routine SLAED8.
00065 *>
00066 *>       The second stage consists of calculating the updated
00067 *>       eigenvalues. This is done by finding the roots of the secular
00068 *>       equation via the routine SLAED4 (as called by SLAED9).
00069 *>       This routine also calculates the eigenvectors of the current
00070 *>       problem.
00071 *>
00072 *>       The final stage consists of computing the updated eigenvectors
00073 *>       directly using the updated eigenvalues.  The eigenvectors for
00074 *>       the current problem are multiplied with the eigenvectors from
00075 *>       the overall problem.
00076 *> \endverbatim
00077 *
00078 *  Arguments:
00079 *  ==========
00080 *
00081 *> \param[in] ICOMPQ
00082 *> \verbatim
00083 *>          ICOMPQ is INTEGER
00084 *>          = 0:  Compute eigenvalues only.
00085 *>          = 1:  Compute eigenvectors of original dense symmetric matrix
00086 *>                also.  On entry, Q contains the orthogonal matrix used
00087 *>                to reduce the original matrix to tridiagonal form.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] N
00091 *> \verbatim
00092 *>          N is INTEGER
00093 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] QSIZ
00097 *> \verbatim
00098 *>          QSIZ is INTEGER
00099 *>         The dimension of the orthogonal matrix used to reduce
00100 *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] TLVLS
00104 *> \verbatim
00105 *>          TLVLS is INTEGER
00106 *>         The total number of merging levels in the overall divide and
00107 *>         conquer tree.
00108 *> \endverbatim
00109 *>
00110 *> \param[in] CURLVL
00111 *> \verbatim
00112 *>          CURLVL is INTEGER
00113 *>         The current level in the overall merge routine,
00114 *>         0 <= CURLVL <= TLVLS.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] CURPBM
00118 *> \verbatim
00119 *>          CURPBM is INTEGER
00120 *>         The current problem in the current level in the overall
00121 *>         merge routine (counting from upper left to lower right).
00122 *> \endverbatim
00123 *>
00124 *> \param[in,out] D
00125 *> \verbatim
00126 *>          D is REAL array, dimension (N)
00127 *>         On entry, the eigenvalues of the rank-1-perturbed matrix.
00128 *>         On exit, the eigenvalues of the repaired matrix.
00129 *> \endverbatim
00130 *>
00131 *> \param[in,out] Q
00132 *> \verbatim
00133 *>          Q is REAL array, dimension (LDQ, N)
00134 *>         On entry, the eigenvectors of the rank-1-perturbed matrix.
00135 *>         On exit, the eigenvectors of the repaired tridiagonal matrix.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] LDQ
00139 *> \verbatim
00140 *>          LDQ is INTEGER
00141 *>         The leading dimension of the array Q.  LDQ >= max(1,N).
00142 *> \endverbatim
00143 *>
00144 *> \param[out] INDXQ
00145 *> \verbatim
00146 *>          INDXQ is INTEGER array, dimension (N)
00147 *>         The permutation which will reintegrate the subproblem just
00148 *>         solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
00149 *>         will be in ascending order.
00150 *> \endverbatim
00151 *>
00152 *> \param[in] RHO
00153 *> \verbatim
00154 *>          RHO is REAL
00155 *>         The subdiagonal element used to create the rank-1
00156 *>         modification.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] CUTPNT
00160 *> \verbatim
00161 *>          CUTPNT is INTEGER
00162 *>         Contains the location of the last eigenvalue in the leading
00163 *>         sub-matrix.  min(1,N) <= CUTPNT <= N.
00164 *> \endverbatim
00165 *>
00166 *> \param[in,out] QSTORE
00167 *> \verbatim
00168 *>          QSTORE is REAL array, dimension (N**2+1)
00169 *>         Stores eigenvectors of submatrices encountered during
00170 *>         divide and conquer, packed together. QPTR points to
00171 *>         beginning of the submatrices.
00172 *> \endverbatim
00173 *>
00174 *> \param[in,out] QPTR
00175 *> \verbatim
00176 *>          QPTR is INTEGER array, dimension (N+2)
00177 *>         List of indices pointing to beginning of submatrices stored
00178 *>         in QSTORE. The submatrices are numbered starting at the
00179 *>         bottom left of the divide and conquer tree, from left to
00180 *>         right and bottom to top.
00181 *> \endverbatim
00182 *>
00183 *> \param[in] PRMPTR
00184 *> \verbatim
00185 *>          PRMPTR is INTEGER array, dimension (N lg N)
00186 *>         Contains a list of pointers which indicate where in PERM a
00187 *>         level's permutation is stored.  PRMPTR(i+1) - PRMPTR(i)
00188 *>         indicates the size of the permutation and also the size of
00189 *>         the full, non-deflated problem.
00190 *> \endverbatim
00191 *>
00192 *> \param[in] PERM
00193 *> \verbatim
00194 *>          PERM is INTEGER array, dimension (N lg N)
00195 *>         Contains the permutations (from deflation and sorting) to be
00196 *>         applied to each eigenblock.
00197 *> \endverbatim
00198 *>
00199 *> \param[in] GIVPTR
00200 *> \verbatim
00201 *>          GIVPTR is INTEGER array, dimension (N lg N)
00202 *>         Contains a list of pointers which indicate where in GIVCOL a
00203 *>         level's Givens rotations are stored.  GIVPTR(i+1) - GIVPTR(i)
00204 *>         indicates the number of Givens rotations.
00205 *> \endverbatim
00206 *>
00207 *> \param[in] GIVCOL
00208 *> \verbatim
00209 *>          GIVCOL is INTEGER array, dimension (2, N lg N)
00210 *>         Each pair of numbers indicates a pair of columns to take place
00211 *>         in a Givens rotation.
00212 *> \endverbatim
00213 *>
00214 *> \param[in] GIVNUM
00215 *> \verbatim
00216 *>          GIVNUM is REAL array, dimension (2, N lg N)
00217 *>         Each number indicates the S value to be used in the
00218 *>         corresponding Givens rotation.
00219 *> \endverbatim
00220 *>
00221 *> \param[out] WORK
00222 *> \verbatim
00223 *>          WORK is REAL array, dimension (3*N+2*QSIZ*N)
00224 *> \endverbatim
00225 *>
00226 *> \param[out] IWORK
00227 *> \verbatim
00228 *>          IWORK is INTEGER array, dimension (4*N)
00229 *> \endverbatim
00230 *>
00231 *> \param[out] INFO
00232 *> \verbatim
00233 *>          INFO is INTEGER
00234 *>          = 0:  successful exit.
00235 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00236 *>          > 0:  if INFO = 1, an eigenvalue did not converge
00237 *> \endverbatim
00238 *
00239 *  Authors:
00240 *  ========
00241 *
00242 *> \author Univ. of Tennessee 
00243 *> \author Univ. of California Berkeley 
00244 *> \author Univ. of Colorado Denver 
00245 *> \author NAG Ltd. 
00246 *
00247 *> \date November 2011
00248 *
00249 *> \ingroup auxOTHERcomputational
00250 *
00251 *> \par Contributors:
00252 *  ==================
00253 *>
00254 *> Jeff Rutter, Computer Science Division, University of California
00255 *> at Berkeley, USA
00256 *
00257 *  =====================================================================
00258       SUBROUTINE SLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
00259      $                   LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
00260      $                   PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
00261      $                   INFO )
00262 *
00263 *  -- LAPACK computational routine (version 3.4.0) --
00264 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00265 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00266 *     November 2011
00267 *
00268 *     .. Scalar Arguments ..
00269       INTEGER            CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
00270      $                   QSIZ, TLVLS
00271       REAL               RHO
00272 *     ..
00273 *     .. Array Arguments ..
00274       INTEGER            GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
00275      $                   IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
00276       REAL               D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
00277      $                   QSTORE( * ), WORK( * )
00278 *     ..
00279 *
00280 *  =====================================================================
00281 *
00282 *     .. Parameters ..
00283       REAL               ONE, ZERO
00284       PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
00285 *     ..
00286 *     .. Local Scalars ..
00287       INTEGER            COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP,
00288      $                   IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR
00289 *     ..
00290 *     .. External Subroutines ..
00291       EXTERNAL           SGEMM, SLAED8, SLAED9, SLAEDA, SLAMRG, XERBLA
00292 *     ..
00293 *     .. Intrinsic Functions ..
00294       INTRINSIC          MAX, MIN
00295 *     ..
00296 *     .. Executable Statements ..
00297 *
00298 *     Test the input parameters.
00299 *
00300       INFO = 0
00301 *
00302       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
00303          INFO = -1
00304       ELSE IF( N.LT.0 ) THEN
00305          INFO = -2
00306       ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
00307          INFO = -4
00308       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00309          INFO = -9
00310       ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
00311          INFO = -12
00312       END IF
00313       IF( INFO.NE.0 ) THEN
00314          CALL XERBLA( 'SLAED7', -INFO )
00315          RETURN
00316       END IF
00317 *
00318 *     Quick return if possible
00319 *
00320       IF( N.EQ.0 )
00321      $   RETURN
00322 *
00323 *     The following values are for bookkeeping purposes only.  They are
00324 *     integer pointers which indicate the portion of the workspace
00325 *     used by a particular array in SLAED8 and SLAED9.
00326 *
00327       IF( ICOMPQ.EQ.1 ) THEN
00328          LDQ2 = QSIZ
00329       ELSE
00330          LDQ2 = N
00331       END IF
00332 *
00333       IZ = 1
00334       IDLMDA = IZ + N
00335       IW = IDLMDA + N
00336       IQ2 = IW + N
00337       IS = IQ2 + N*LDQ2
00338 *
00339       INDX = 1
00340       INDXC = INDX + N
00341       COLTYP = INDXC + N
00342       INDXP = COLTYP + N
00343 *
00344 *     Form the z-vector which consists of the last row of Q_1 and the
00345 *     first row of Q_2.
00346 *
00347       PTR = 1 + 2**TLVLS
00348       DO 10 I = 1, CURLVL - 1
00349          PTR = PTR + 2**( TLVLS-I )
00350    10 CONTINUE
00351       CURR = PTR + CURPBM
00352       CALL SLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
00353      $             GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ),
00354      $             WORK( IZ+N ), INFO )
00355 *
00356 *     When solving the final problem, we no longer need the stored data,
00357 *     so we will overwrite the data from this level onto the previously
00358 *     used storage space.
00359 *
00360       IF( CURLVL.EQ.TLVLS ) THEN
00361          QPTR( CURR ) = 1
00362          PRMPTR( CURR ) = 1
00363          GIVPTR( CURR ) = 1
00364       END IF
00365 *
00366 *     Sort and Deflate eigenvalues.
00367 *
00368       CALL SLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT,
00369      $             WORK( IZ ), WORK( IDLMDA ), WORK( IQ2 ), LDQ2,
00370      $             WORK( IW ), PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
00371      $             GIVCOL( 1, GIVPTR( CURR ) ),
00372      $             GIVNUM( 1, GIVPTR( CURR ) ), IWORK( INDXP ),
00373      $             IWORK( INDX ), INFO )
00374       PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
00375       GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
00376 *
00377 *     Solve Secular Equation.
00378 *
00379       IF( K.NE.0 ) THEN
00380          CALL SLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ),
00381      $                WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO )
00382          IF( INFO.NE.0 )
00383      $      GO TO 30
00384          IF( ICOMPQ.EQ.1 ) THEN
00385             CALL SGEMM( 'N', 'N', QSIZ, K, K, ONE, WORK( IQ2 ), LDQ2,
00386      $                  QSTORE( QPTR( CURR ) ), K, ZERO, Q, LDQ )
00387          END IF
00388          QPTR( CURR+1 ) = QPTR( CURR ) + K**2
00389 *
00390 *     Prepare the INDXQ sorting permutation.
00391 *
00392          N1 = K
00393          N2 = N - K
00394          CALL SLAMRG( N1, N2, D, 1, -1, INDXQ )
00395       ELSE
00396          QPTR( CURR+1 ) = QPTR( CURR )
00397          DO 20 I = 1, N
00398             INDXQ( I ) = I
00399    20    CONTINUE
00400       END IF
00401 *
00402    30 CONTINUE
00403       RETURN
00404 *
00405 *     End of SLAED7
00406 *
00407       END
 All Files Functions