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