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