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