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