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