![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLAED8 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLAED8 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claed8.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claed8.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claed8.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, 00022 * Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, 00023 * GIVCOL, GIVNUM, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ 00027 * REAL RHO 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), 00031 * $ INDXQ( * ), PERM( * ) 00032 * REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), 00033 * $ Z( * ) 00034 * COMPLEX Q( LDQ, * ), Q2( LDQ2, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> CLAED8 merges the two sets of eigenvalues together into a single 00044 *> sorted set. Then it tries to deflate the size of the problem. 00045 *> There are two ways in which deflation can occur: when two or more 00046 *> eigenvalues are close together or if there is a tiny element in the 00047 *> Z vector. For each such occurrence the order of the related secular 00048 *> equation problem is reduced by one. 00049 *> \endverbatim 00050 * 00051 * Arguments: 00052 * ========== 00053 * 00054 *> \param[out] K 00055 *> \verbatim 00056 *> K is INTEGER 00057 *> Contains the number of non-deflated eigenvalues. 00058 *> This is the order of the related secular equation. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] N 00062 *> \verbatim 00063 *> N is INTEGER 00064 *> The dimension of the symmetric tridiagonal matrix. N >= 0. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] QSIZ 00068 *> \verbatim 00069 *> QSIZ is INTEGER 00070 *> The dimension of the unitary matrix used to reduce 00071 *> the dense or band matrix to tridiagonal form. 00072 *> QSIZ >= N if ICOMPQ = 1. 00073 *> \endverbatim 00074 *> 00075 *> \param[in,out] Q 00076 *> \verbatim 00077 *> Q is COMPLEX array, dimension (LDQ,N) 00078 *> On entry, Q contains the eigenvectors of the partially solved 00079 *> system which has been previously updated in matrix 00080 *> multiplies with other partially solved eigensystems. 00081 *> On exit, Q contains the trailing (N-K) updated eigenvectors 00082 *> (those which were deflated) in its last N-K columns. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] LDQ 00086 *> \verbatim 00087 *> LDQ is INTEGER 00088 *> The leading dimension of the array Q. LDQ >= max( 1, N ). 00089 *> \endverbatim 00090 *> 00091 *> \param[in,out] D 00092 *> \verbatim 00093 *> D is REAL array, dimension (N) 00094 *> On entry, D contains the eigenvalues of the two submatrices to 00095 *> be combined. On exit, D contains the trailing (N-K) updated 00096 *> eigenvalues (those which were deflated) sorted into increasing 00097 *> order. 00098 *> \endverbatim 00099 *> 00100 *> \param[in,out] RHO 00101 *> \verbatim 00102 *> RHO is REAL 00103 *> Contains the off diagonal element associated with the rank-1 00104 *> cut which originally split the two submatrices which are now 00105 *> being recombined. RHO is modified during the computation to 00106 *> the value required by SLAED3. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] CUTPNT 00110 *> \verbatim 00111 *> CUTPNT is INTEGER 00112 *> Contains the location of the last eigenvalue in the leading 00113 *> sub-matrix. MIN(1,N) <= CUTPNT <= N. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] Z 00117 *> \verbatim 00118 *> Z is REAL array, dimension (N) 00119 *> On input this vector contains the updating vector (the last 00120 *> row of the first sub-eigenvector matrix and the first row of 00121 *> the second sub-eigenvector matrix). The contents of Z are 00122 *> destroyed during the updating process. 00123 *> \endverbatim 00124 *> 00125 *> \param[out] DLAMDA 00126 *> \verbatim 00127 *> DLAMDA is REAL array, dimension (N) 00128 *> Contains a copy of the first K eigenvalues which will be used 00129 *> by SLAED3 to form the secular equation. 00130 *> \endverbatim 00131 *> 00132 *> \param[out] Q2 00133 *> \verbatim 00134 *> Q2 is COMPLEX array, dimension (LDQ2,N) 00135 *> If ICOMPQ = 0, Q2 is not referenced. Otherwise, 00136 *> Contains a copy of the first K eigenvectors which will be used 00137 *> by SLAED7 in a matrix multiply (SGEMM) to update the new 00138 *> eigenvectors. 00139 *> \endverbatim 00140 *> 00141 *> \param[in] LDQ2 00142 *> \verbatim 00143 *> LDQ2 is INTEGER 00144 *> The leading dimension of the array Q2. LDQ2 >= max( 1, N ). 00145 *> \endverbatim 00146 *> 00147 *> \param[out] W 00148 *> \verbatim 00149 *> W is REAL array, dimension (N) 00150 *> This will hold the first k values of the final 00151 *> deflation-altered z-vector and will be passed to SLAED3. 00152 *> \endverbatim 00153 *> 00154 *> \param[out] INDXP 00155 *> \verbatim 00156 *> INDXP is INTEGER array, dimension (N) 00157 *> This will contain the permutation used to place deflated 00158 *> values of D at the end of the array. On output INDXP(1:K) 00159 *> points to the nondeflated D-values and INDXP(K+1:N) 00160 *> points to the deflated eigenvalues. 00161 *> \endverbatim 00162 *> 00163 *> \param[out] INDX 00164 *> \verbatim 00165 *> INDX is INTEGER array, dimension (N) 00166 *> This will contain the permutation used to sort the contents of 00167 *> D into ascending order. 00168 *> \endverbatim 00169 *> 00170 *> \param[in] INDXQ 00171 *> \verbatim 00172 *> INDXQ is INTEGER array, dimension (N) 00173 *> This contains the permutation which separately sorts the two 00174 *> sub-problems in D into ascending order. Note that elements in 00175 *> the second half of this permutation must first have CUTPNT 00176 *> added to their values in order to be accurate. 00177 *> \endverbatim 00178 *> 00179 *> \param[out] PERM 00180 *> \verbatim 00181 *> PERM is INTEGER array, dimension (N) 00182 *> Contains the permutations (from deflation and sorting) to be 00183 *> applied to each eigenblock. 00184 *> \endverbatim 00185 *> 00186 *> \param[out] GIVPTR 00187 *> \verbatim 00188 *> GIVPTR is INTEGER 00189 *> Contains the number of Givens rotations which took place in 00190 *> this subproblem. 00191 *> \endverbatim 00192 *> 00193 *> \param[out] GIVCOL 00194 *> \verbatim 00195 *> GIVCOL is INTEGER array, dimension (2, N) 00196 *> Each pair of numbers indicates a pair of columns to take place 00197 *> in a Givens rotation. 00198 *> \endverbatim 00199 *> 00200 *> \param[out] GIVNUM 00201 *> \verbatim 00202 *> GIVNUM is REAL array, dimension (2, N) 00203 *> Each number indicates the S value to be used in the 00204 *> corresponding Givens rotation. 00205 *> \endverbatim 00206 *> 00207 *> \param[out] INFO 00208 *> \verbatim 00209 *> INFO is INTEGER 00210 *> = 0: successful exit. 00211 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00212 *> \endverbatim 00213 * 00214 * Authors: 00215 * ======== 00216 * 00217 *> \author Univ. of Tennessee 00218 *> \author Univ. of California Berkeley 00219 *> \author Univ. of Colorado Denver 00220 *> \author NAG Ltd. 00221 * 00222 *> \date November 2011 00223 * 00224 *> \ingroup complexOTHERcomputational 00225 * 00226 * ===================================================================== 00227 SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, 00228 $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, 00229 $ GIVCOL, GIVNUM, INFO ) 00230 * 00231 * -- LAPACK computational routine (version 3.4.0) -- 00232 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00233 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00234 * November 2011 00235 * 00236 * .. Scalar Arguments .. 00237 INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ 00238 REAL RHO 00239 * .. 00240 * .. Array Arguments .. 00241 INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), 00242 $ INDXQ( * ), PERM( * ) 00243 REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), 00244 $ Z( * ) 00245 COMPLEX Q( LDQ, * ), Q2( LDQ2, * ) 00246 * .. 00247 * 00248 * ===================================================================== 00249 * 00250 * .. Parameters .. 00251 REAL MONE, ZERO, ONE, TWO, EIGHT 00252 PARAMETER ( MONE = -1.0E0, ZERO = 0.0E0, ONE = 1.0E0, 00253 $ TWO = 2.0E0, EIGHT = 8.0E0 ) 00254 * .. 00255 * .. Local Scalars .. 00256 INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2 00257 REAL C, EPS, S, T, TAU, TOL 00258 * .. 00259 * .. External Functions .. 00260 INTEGER ISAMAX 00261 REAL SLAMCH, SLAPY2 00262 EXTERNAL ISAMAX, SLAMCH, SLAPY2 00263 * .. 00264 * .. External Subroutines .. 00265 EXTERNAL CCOPY, CLACPY, CSROT, SCOPY, SLAMRG, SSCAL, 00266 $ XERBLA 00267 * .. 00268 * .. Intrinsic Functions .. 00269 INTRINSIC ABS, MAX, MIN, SQRT 00270 * .. 00271 * .. Executable Statements .. 00272 * 00273 * Test the input parameters. 00274 * 00275 INFO = 0 00276 * 00277 IF( N.LT.0 ) THEN 00278 INFO = -2 00279 ELSE IF( QSIZ.LT.N ) THEN 00280 INFO = -3 00281 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00282 INFO = -5 00283 ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN 00284 INFO = -8 00285 ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN 00286 INFO = -12 00287 END IF 00288 IF( INFO.NE.0 ) THEN 00289 CALL XERBLA( 'CLAED8', -INFO ) 00290 RETURN 00291 END IF 00292 * 00293 * Need to initialize GIVPTR to O here in case of quick exit 00294 * to prevent an unspecified code behavior (usually sigfault) 00295 * when IWORK array on entry to *stedc is not zeroed 00296 * (or at least some IWORK entries which used in *laed7 for GIVPTR). 00297 * 00298 GIVPTR = 0 00299 * 00300 * Quick return if possible 00301 * 00302 IF( N.EQ.0 ) 00303 $ RETURN 00304 * 00305 N1 = CUTPNT 00306 N2 = N - N1 00307 N1P1 = N1 + 1 00308 * 00309 IF( RHO.LT.ZERO ) THEN 00310 CALL SSCAL( N2, MONE, Z( N1P1 ), 1 ) 00311 END IF 00312 * 00313 * Normalize z so that norm(z) = 1 00314 * 00315 T = ONE / SQRT( TWO ) 00316 DO 10 J = 1, N 00317 INDX( J ) = J 00318 10 CONTINUE 00319 CALL SSCAL( N, T, Z, 1 ) 00320 RHO = ABS( TWO*RHO ) 00321 * 00322 * Sort the eigenvalues into increasing order 00323 * 00324 DO 20 I = CUTPNT + 1, N 00325 INDXQ( I ) = INDXQ( I ) + CUTPNT 00326 20 CONTINUE 00327 DO 30 I = 1, N 00328 DLAMDA( I ) = D( INDXQ( I ) ) 00329 W( I ) = Z( INDXQ( I ) ) 00330 30 CONTINUE 00331 I = 1 00332 J = CUTPNT + 1 00333 CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) 00334 DO 40 I = 1, N 00335 D( I ) = DLAMDA( INDX( I ) ) 00336 Z( I ) = W( INDX( I ) ) 00337 40 CONTINUE 00338 * 00339 * Calculate the allowable deflation tolerance 00340 * 00341 IMAX = ISAMAX( N, Z, 1 ) 00342 JMAX = ISAMAX( N, D, 1 ) 00343 EPS = SLAMCH( 'Epsilon' ) 00344 TOL = EIGHT*EPS*ABS( D( JMAX ) ) 00345 * 00346 * If the rank-1 modifier is small enough, no more needs to be done 00347 * -- except to reorganize Q so that its columns correspond with the 00348 * elements in D. 00349 * 00350 IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN 00351 K = 0 00352 DO 50 J = 1, N 00353 PERM( J ) = INDXQ( INDX( J ) ) 00354 CALL CCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 00355 50 CONTINUE 00356 CALL CLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), LDQ ) 00357 RETURN 00358 END IF 00359 * 00360 * If there are multiple eigenvalues then the problem deflates. Here 00361 * the number of equal eigenvalues are found. As each equal 00362 * eigenvalue is found, an elementary reflector is computed to rotate 00363 * the corresponding eigensubspace so that the corresponding 00364 * components of Z are zero in this new basis. 00365 * 00366 K = 0 00367 K2 = N + 1 00368 DO 60 J = 1, N 00369 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 00370 * 00371 * Deflate due to small z component. 00372 * 00373 K2 = K2 - 1 00374 INDXP( K2 ) = J 00375 IF( J.EQ.N ) 00376 $ GO TO 100 00377 ELSE 00378 JLAM = J 00379 GO TO 70 00380 END IF 00381 60 CONTINUE 00382 70 CONTINUE 00383 J = J + 1 00384 IF( J.GT.N ) 00385 $ GO TO 90 00386 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN 00387 * 00388 * Deflate due to small z component. 00389 * 00390 K2 = K2 - 1 00391 INDXP( K2 ) = J 00392 ELSE 00393 * 00394 * Check if eigenvalues are close enough to allow deflation. 00395 * 00396 S = Z( JLAM ) 00397 C = Z( J ) 00398 * 00399 * Find sqrt(a**2+b**2) without overflow or 00400 * destructive underflow. 00401 * 00402 TAU = SLAPY2( C, S ) 00403 T = D( J ) - D( JLAM ) 00404 C = C / TAU 00405 S = -S / TAU 00406 IF( ABS( T*C*S ).LE.TOL ) THEN 00407 * 00408 * Deflation is possible. 00409 * 00410 Z( J ) = TAU 00411 Z( JLAM ) = ZERO 00412 * 00413 * Record the appropriate Givens rotation 00414 * 00415 GIVPTR = GIVPTR + 1 00416 GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) ) 00417 GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) ) 00418 GIVNUM( 1, GIVPTR ) = C 00419 GIVNUM( 2, GIVPTR ) = S 00420 CALL CSROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1, 00421 $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S ) 00422 T = D( JLAM )*C*C + D( J )*S*S 00423 D( J ) = D( JLAM )*S*S + D( J )*C*C 00424 D( JLAM ) = T 00425 K2 = K2 - 1 00426 I = 1 00427 80 CONTINUE 00428 IF( K2+I.LE.N ) THEN 00429 IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN 00430 INDXP( K2+I-1 ) = INDXP( K2+I ) 00431 INDXP( K2+I ) = JLAM 00432 I = I + 1 00433 GO TO 80 00434 ELSE 00435 INDXP( K2+I-1 ) = JLAM 00436 END IF 00437 ELSE 00438 INDXP( K2+I-1 ) = JLAM 00439 END IF 00440 JLAM = J 00441 ELSE 00442 K = K + 1 00443 W( K ) = Z( JLAM ) 00444 DLAMDA( K ) = D( JLAM ) 00445 INDXP( K ) = JLAM 00446 JLAM = J 00447 END IF 00448 END IF 00449 GO TO 70 00450 90 CONTINUE 00451 * 00452 * Record the last eigenvalue. 00453 * 00454 K = K + 1 00455 W( K ) = Z( JLAM ) 00456 DLAMDA( K ) = D( JLAM ) 00457 INDXP( K ) = JLAM 00458 * 00459 100 CONTINUE 00460 * 00461 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA 00462 * and Q2 respectively. The eigenvalues/vectors which were not 00463 * deflated go into the first K slots of DLAMDA and Q2 respectively, 00464 * while those which were deflated go into the last N - K slots. 00465 * 00466 DO 110 J = 1, N 00467 JP = INDXP( J ) 00468 DLAMDA( J ) = D( JP ) 00469 PERM( J ) = INDXQ( INDX( JP ) ) 00470 CALL CCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) 00471 110 CONTINUE 00472 * 00473 * The deflated eigenvalues and their corresponding vectors go back 00474 * into the last N - K slots of D and Q respectively. 00475 * 00476 IF( K.LT.N ) THEN 00477 CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) 00478 CALL CLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ), 00479 $ LDQ ) 00480 END IF 00481 * 00482 RETURN 00483 * 00484 * End of CLAED8 00485 * 00486 END