![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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