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