1 SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
2 $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
3 $ GIVCOL, GIVNUM, INFO )
4 *
5 * -- LAPACK routine (version 3.2.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * June 2010
9 *
10 * .. Scalar Arguments ..
11 INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
12 DOUBLE PRECISION RHO
13 * ..
14 * .. Array Arguments ..
15 INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
16 $ INDXQ( * ), PERM( * )
17 DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
18 $ Z( * )
19 COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZLAED8 merges the two sets of eigenvalues together into a single
26 * sorted set. Then it tries to deflate the size of the problem.
27 * There are two ways in which deflation can occur: when two or more
28 * eigenvalues are close together or if there is a tiny element in the
29 * Z vector. For each such occurrence the order of the related secular
30 * equation problem is reduced by one.
31 *
32 * Arguments
33 * =========
34 *
35 * K (output) INTEGER
36 * Contains the number of non-deflated eigenvalues.
37 * This is the order of the related secular equation.
38 *
39 * N (input) INTEGER
40 * The dimension of the symmetric tridiagonal matrix. N >= 0.
41 *
42 * QSIZ (input) INTEGER
43 * The dimension of the unitary matrix used to reduce
44 * the dense or band matrix to tridiagonal form.
45 * QSIZ >= N if ICOMPQ = 1.
46 *
47 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
48 * On entry, Q contains the eigenvectors of the partially solved
49 * system which has been previously updated in matrix
50 * multiplies with other partially solved eigensystems.
51 * On exit, Q contains the trailing (N-K) updated eigenvectors
52 * (those which were deflated) in its last N-K columns.
53 *
54 * LDQ (input) INTEGER
55 * The leading dimension of the array Q. LDQ >= max( 1, N ).
56 *
57 * D (input/output) DOUBLE PRECISION array, dimension (N)
58 * On entry, D contains the eigenvalues of the two submatrices to
59 * be combined. On exit, D contains the trailing (N-K) updated
60 * eigenvalues (those which were deflated) sorted into increasing
61 * order.
62 *
63 * RHO (input/output) DOUBLE PRECISION
64 * Contains the off diagonal element associated with the rank-1
65 * cut which originally split the two submatrices which are now
66 * being recombined. RHO is modified during the computation to
67 * the value required by DLAED3.
68 *
69 * CUTPNT (input) INTEGER
70 * Contains the location of the last eigenvalue in the leading
71 * sub-matrix. MIN(1,N) <= CUTPNT <= N.
72 *
73 * Z (input) DOUBLE PRECISION array, dimension (N)
74 * On input this vector contains the updating vector (the last
75 * row of the first sub-eigenvector matrix and the first row of
76 * the second sub-eigenvector matrix). The contents of Z are
77 * destroyed during the updating process.
78 *
79 * DLAMDA (output) DOUBLE PRECISION array, dimension (N)
80 * Contains a copy of the first K eigenvalues which will be used
81 * by DLAED3 to form the secular equation.
82 *
83 * Q2 (output) COMPLEX*16 array, dimension (LDQ2,N)
84 * If ICOMPQ = 0, Q2 is not referenced. Otherwise,
85 * Contains a copy of the first K eigenvectors which will be used
86 * by DLAED7 in a matrix multiply (DGEMM) to update the new
87 * eigenvectors.
88 *
89 * LDQ2 (input) INTEGER
90 * The leading dimension of the array Q2. LDQ2 >= max( 1, N ).
91 *
92 * W (output) DOUBLE PRECISION array, dimension (N)
93 * This will hold the first k values of the final
94 * deflation-altered z-vector and will be passed to DLAED3.
95 *
96 * INDXP (workspace) INTEGER array, dimension (N)
97 * This will contain the permutation used to place deflated
98 * values of D at the end of the array. On output INDXP(1:K)
99 * points to the nondeflated D-values and INDXP(K+1:N)
100 * points to the deflated eigenvalues.
101 *
102 * INDX (workspace) INTEGER array, dimension (N)
103 * This will contain the permutation used to sort the contents of
104 * D into ascending order.
105 *
106 * INDXQ (input) INTEGER array, dimension (N)
107 * This contains the permutation which separately sorts the two
108 * sub-problems in D into ascending order. Note that elements in
109 * the second half of this permutation must first have CUTPNT
110 * added to their values in order to be accurate.
111 *
112 * PERM (output) INTEGER array, dimension (N)
113 * Contains the permutations (from deflation and sorting) to be
114 * applied to each eigenblock.
115 *
116 * GIVPTR (output) INTEGER
117 * Contains the number of Givens rotations which took place in
118 * this subproblem.
119 *
120 * GIVCOL (output) INTEGER array, dimension (2, N)
121 * Each pair of numbers indicates a pair of columns to take place
122 * in a Givens rotation.
123 *
124 * GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
125 * Each number indicates the S value to be used in the
126 * corresponding Givens rotation.
127 *
128 * INFO (output) INTEGER
129 * = 0: successful exit.
130 * < 0: if INFO = -i, the i-th argument had an illegal value.
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135 DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
136 PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
137 $ TWO = 2.0D0, EIGHT = 8.0D0 )
138 * ..
139 * .. Local Scalars ..
140 INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
141 DOUBLE PRECISION C, EPS, S, T, TAU, TOL
142 * ..
143 * .. External Functions ..
144 INTEGER IDAMAX
145 DOUBLE PRECISION DLAMCH, DLAPY2
146 EXTERNAL IDAMAX, DLAMCH, DLAPY2
147 * ..
148 * .. External Subroutines ..
149 EXTERNAL DCOPY, DLAMRG, DSCAL, XERBLA, ZCOPY, ZDROT,
150 $ ZLACPY
151 * ..
152 * .. Intrinsic Functions ..
153 INTRINSIC ABS, MAX, MIN, SQRT
154 * ..
155 * .. Executable Statements ..
156 *
157 * Test the input parameters.
158 *
159 INFO = 0
160 *
161 IF( N.LT.0 ) THEN
162 INFO = -2
163 ELSE IF( QSIZ.LT.N ) THEN
164 INFO = -3
165 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
166 INFO = -5
167 ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
168 INFO = -8
169 ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
170 INFO = -12
171 END IF
172 IF( INFO.NE.0 ) THEN
173 CALL XERBLA( 'ZLAED8', -INFO )
174 RETURN
175 END IF
176 *
177 * Need to initialize GIVPTR to O here in case of quick exit
178 * to prevent an unspecified code behavior (usually sigfault)
179 * when IWORK array on entry to *stedc is not zeroed
180 * (or at least some IWORK entries which used in *laed7 for GIVPTR).
181 *
182 GIVPTR = 0
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 N1 = CUTPNT
190 N2 = N - N1
191 N1P1 = N1 + 1
192 *
193 IF( RHO.LT.ZERO ) THEN
194 CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
195 END IF
196 *
197 * Normalize z so that norm(z) = 1
198 *
199 T = ONE / SQRT( TWO )
200 DO 10 J = 1, N
201 INDX( J ) = J
202 10 CONTINUE
203 CALL DSCAL( N, T, Z, 1 )
204 RHO = ABS( TWO*RHO )
205 *
206 * Sort the eigenvalues into increasing order
207 *
208 DO 20 I = CUTPNT + 1, N
209 INDXQ( I ) = INDXQ( I ) + CUTPNT
210 20 CONTINUE
211 DO 30 I = 1, N
212 DLAMDA( I ) = D( INDXQ( I ) )
213 W( I ) = Z( INDXQ( I ) )
214 30 CONTINUE
215 I = 1
216 J = CUTPNT + 1
217 CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
218 DO 40 I = 1, N
219 D( I ) = DLAMDA( INDX( I ) )
220 Z( I ) = W( INDX( I ) )
221 40 CONTINUE
222 *
223 * Calculate the allowable deflation tolerance
224 *
225 IMAX = IDAMAX( N, Z, 1 )
226 JMAX = IDAMAX( N, D, 1 )
227 EPS = DLAMCH( 'Epsilon' )
228 TOL = EIGHT*EPS*ABS( D( JMAX ) )
229 *
230 * If the rank-1 modifier is small enough, no more needs to be done
231 * -- except to reorganize Q so that its columns correspond with the
232 * elements in D.
233 *
234 IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
235 K = 0
236 DO 50 J = 1, N
237 PERM( J ) = INDXQ( INDX( J ) )
238 CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
239 50 CONTINUE
240 CALL ZLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), LDQ )
241 RETURN
242 END IF
243 *
244 * If there are multiple eigenvalues then the problem deflates. Here
245 * the number of equal eigenvalues are found. As each equal
246 * eigenvalue is found, an elementary reflector is computed to rotate
247 * the corresponding eigensubspace so that the corresponding
248 * components of Z are zero in this new basis.
249 *
250 K = 0
251 K2 = N + 1
252 DO 60 J = 1, N
253 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
254 *
255 * Deflate due to small z component.
256 *
257 K2 = K2 - 1
258 INDXP( K2 ) = J
259 IF( J.EQ.N )
260 $ GO TO 100
261 ELSE
262 JLAM = J
263 GO TO 70
264 END IF
265 60 CONTINUE
266 70 CONTINUE
267 J = J + 1
268 IF( J.GT.N )
269 $ GO TO 90
270 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
271 *
272 * Deflate due to small z component.
273 *
274 K2 = K2 - 1
275 INDXP( K2 ) = J
276 ELSE
277 *
278 * Check if eigenvalues are close enough to allow deflation.
279 *
280 S = Z( JLAM )
281 C = Z( J )
282 *
283 * Find sqrt(a**2+b**2) without overflow or
284 * destructive underflow.
285 *
286 TAU = DLAPY2( C, S )
287 T = D( J ) - D( JLAM )
288 C = C / TAU
289 S = -S / TAU
290 IF( ABS( T*C*S ).LE.TOL ) THEN
291 *
292 * Deflation is possible.
293 *
294 Z( J ) = TAU
295 Z( JLAM ) = ZERO
296 *
297 * Record the appropriate Givens rotation
298 *
299 GIVPTR = GIVPTR + 1
300 GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
301 GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
302 GIVNUM( 1, GIVPTR ) = C
303 GIVNUM( 2, GIVPTR ) = S
304 CALL ZDROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
305 $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
306 T = D( JLAM )*C*C + D( J )*S*S
307 D( J ) = D( JLAM )*S*S + D( J )*C*C
308 D( JLAM ) = T
309 K2 = K2 - 1
310 I = 1
311 80 CONTINUE
312 IF( K2+I.LE.N ) THEN
313 IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
314 INDXP( K2+I-1 ) = INDXP( K2+I )
315 INDXP( K2+I ) = JLAM
316 I = I + 1
317 GO TO 80
318 ELSE
319 INDXP( K2+I-1 ) = JLAM
320 END IF
321 ELSE
322 INDXP( K2+I-1 ) = JLAM
323 END IF
324 JLAM = J
325 ELSE
326 K = K + 1
327 W( K ) = Z( JLAM )
328 DLAMDA( K ) = D( JLAM )
329 INDXP( K ) = JLAM
330 JLAM = J
331 END IF
332 END IF
333 GO TO 70
334 90 CONTINUE
335 *
336 * Record the last eigenvalue.
337 *
338 K = K + 1
339 W( K ) = Z( JLAM )
340 DLAMDA( K ) = D( JLAM )
341 INDXP( K ) = JLAM
342 *
343 100 CONTINUE
344 *
345 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
346 * and Q2 respectively. The eigenvalues/vectors which were not
347 * deflated go into the first K slots of DLAMDA and Q2 respectively,
348 * while those which were deflated go into the last N - K slots.
349 *
350 DO 110 J = 1, N
351 JP = INDXP( J )
352 DLAMDA( J ) = D( JP )
353 PERM( J ) = INDXQ( INDX( JP ) )
354 CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
355 110 CONTINUE
356 *
357 * The deflated eigenvalues and their corresponding vectors go back
358 * into the last N - K slots of D and Q respectively.
359 *
360 IF( K.LT.N ) THEN
361 CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
362 CALL ZLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ),
363 $ LDQ )
364 END IF
365 *
366 RETURN
367 *
368 * End of ZLAED8
369 *
370 END
2 $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
3 $ GIVCOL, GIVNUM, INFO )
4 *
5 * -- LAPACK routine (version 3.2.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * June 2010
9 *
10 * .. Scalar Arguments ..
11 INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
12 DOUBLE PRECISION RHO
13 * ..
14 * .. Array Arguments ..
15 INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
16 $ INDXQ( * ), PERM( * )
17 DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
18 $ Z( * )
19 COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * ZLAED8 merges the two sets of eigenvalues together into a single
26 * sorted set. Then it tries to deflate the size of the problem.
27 * There are two ways in which deflation can occur: when two or more
28 * eigenvalues are close together or if there is a tiny element in the
29 * Z vector. For each such occurrence the order of the related secular
30 * equation problem is reduced by one.
31 *
32 * Arguments
33 * =========
34 *
35 * K (output) INTEGER
36 * Contains the number of non-deflated eigenvalues.
37 * This is the order of the related secular equation.
38 *
39 * N (input) INTEGER
40 * The dimension of the symmetric tridiagonal matrix. N >= 0.
41 *
42 * QSIZ (input) INTEGER
43 * The dimension of the unitary matrix used to reduce
44 * the dense or band matrix to tridiagonal form.
45 * QSIZ >= N if ICOMPQ = 1.
46 *
47 * Q (input/output) COMPLEX*16 array, dimension (LDQ,N)
48 * On entry, Q contains the eigenvectors of the partially solved
49 * system which has been previously updated in matrix
50 * multiplies with other partially solved eigensystems.
51 * On exit, Q contains the trailing (N-K) updated eigenvectors
52 * (those which were deflated) in its last N-K columns.
53 *
54 * LDQ (input) INTEGER
55 * The leading dimension of the array Q. LDQ >= max( 1, N ).
56 *
57 * D (input/output) DOUBLE PRECISION array, dimension (N)
58 * On entry, D contains the eigenvalues of the two submatrices to
59 * be combined. On exit, D contains the trailing (N-K) updated
60 * eigenvalues (those which were deflated) sorted into increasing
61 * order.
62 *
63 * RHO (input/output) DOUBLE PRECISION
64 * Contains the off diagonal element associated with the rank-1
65 * cut which originally split the two submatrices which are now
66 * being recombined. RHO is modified during the computation to
67 * the value required by DLAED3.
68 *
69 * CUTPNT (input) INTEGER
70 * Contains the location of the last eigenvalue in the leading
71 * sub-matrix. MIN(1,N) <= CUTPNT <= N.
72 *
73 * Z (input) DOUBLE PRECISION array, dimension (N)
74 * On input this vector contains the updating vector (the last
75 * row of the first sub-eigenvector matrix and the first row of
76 * the second sub-eigenvector matrix). The contents of Z are
77 * destroyed during the updating process.
78 *
79 * DLAMDA (output) DOUBLE PRECISION array, dimension (N)
80 * Contains a copy of the first K eigenvalues which will be used
81 * by DLAED3 to form the secular equation.
82 *
83 * Q2 (output) COMPLEX*16 array, dimension (LDQ2,N)
84 * If ICOMPQ = 0, Q2 is not referenced. Otherwise,
85 * Contains a copy of the first K eigenvectors which will be used
86 * by DLAED7 in a matrix multiply (DGEMM) to update the new
87 * eigenvectors.
88 *
89 * LDQ2 (input) INTEGER
90 * The leading dimension of the array Q2. LDQ2 >= max( 1, N ).
91 *
92 * W (output) DOUBLE PRECISION array, dimension (N)
93 * This will hold the first k values of the final
94 * deflation-altered z-vector and will be passed to DLAED3.
95 *
96 * INDXP (workspace) INTEGER array, dimension (N)
97 * This will contain the permutation used to place deflated
98 * values of D at the end of the array. On output INDXP(1:K)
99 * points to the nondeflated D-values and INDXP(K+1:N)
100 * points to the deflated eigenvalues.
101 *
102 * INDX (workspace) INTEGER array, dimension (N)
103 * This will contain the permutation used to sort the contents of
104 * D into ascending order.
105 *
106 * INDXQ (input) INTEGER array, dimension (N)
107 * This contains the permutation which separately sorts the two
108 * sub-problems in D into ascending order. Note that elements in
109 * the second half of this permutation must first have CUTPNT
110 * added to their values in order to be accurate.
111 *
112 * PERM (output) INTEGER array, dimension (N)
113 * Contains the permutations (from deflation and sorting) to be
114 * applied to each eigenblock.
115 *
116 * GIVPTR (output) INTEGER
117 * Contains the number of Givens rotations which took place in
118 * this subproblem.
119 *
120 * GIVCOL (output) INTEGER array, dimension (2, N)
121 * Each pair of numbers indicates a pair of columns to take place
122 * in a Givens rotation.
123 *
124 * GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
125 * Each number indicates the S value to be used in the
126 * corresponding Givens rotation.
127 *
128 * INFO (output) INTEGER
129 * = 0: successful exit.
130 * < 0: if INFO = -i, the i-th argument had an illegal value.
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135 DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
136 PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
137 $ TWO = 2.0D0, EIGHT = 8.0D0 )
138 * ..
139 * .. Local Scalars ..
140 INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
141 DOUBLE PRECISION C, EPS, S, T, TAU, TOL
142 * ..
143 * .. External Functions ..
144 INTEGER IDAMAX
145 DOUBLE PRECISION DLAMCH, DLAPY2
146 EXTERNAL IDAMAX, DLAMCH, DLAPY2
147 * ..
148 * .. External Subroutines ..
149 EXTERNAL DCOPY, DLAMRG, DSCAL, XERBLA, ZCOPY, ZDROT,
150 $ ZLACPY
151 * ..
152 * .. Intrinsic Functions ..
153 INTRINSIC ABS, MAX, MIN, SQRT
154 * ..
155 * .. Executable Statements ..
156 *
157 * Test the input parameters.
158 *
159 INFO = 0
160 *
161 IF( N.LT.0 ) THEN
162 INFO = -2
163 ELSE IF( QSIZ.LT.N ) THEN
164 INFO = -3
165 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
166 INFO = -5
167 ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
168 INFO = -8
169 ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
170 INFO = -12
171 END IF
172 IF( INFO.NE.0 ) THEN
173 CALL XERBLA( 'ZLAED8', -INFO )
174 RETURN
175 END IF
176 *
177 * Need to initialize GIVPTR to O here in case of quick exit
178 * to prevent an unspecified code behavior (usually sigfault)
179 * when IWORK array on entry to *stedc is not zeroed
180 * (or at least some IWORK entries which used in *laed7 for GIVPTR).
181 *
182 GIVPTR = 0
183 *
184 * Quick return if possible
185 *
186 IF( N.EQ.0 )
187 $ RETURN
188 *
189 N1 = CUTPNT
190 N2 = N - N1
191 N1P1 = N1 + 1
192 *
193 IF( RHO.LT.ZERO ) THEN
194 CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
195 END IF
196 *
197 * Normalize z so that norm(z) = 1
198 *
199 T = ONE / SQRT( TWO )
200 DO 10 J = 1, N
201 INDX( J ) = J
202 10 CONTINUE
203 CALL DSCAL( N, T, Z, 1 )
204 RHO = ABS( TWO*RHO )
205 *
206 * Sort the eigenvalues into increasing order
207 *
208 DO 20 I = CUTPNT + 1, N
209 INDXQ( I ) = INDXQ( I ) + CUTPNT
210 20 CONTINUE
211 DO 30 I = 1, N
212 DLAMDA( I ) = D( INDXQ( I ) )
213 W( I ) = Z( INDXQ( I ) )
214 30 CONTINUE
215 I = 1
216 J = CUTPNT + 1
217 CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
218 DO 40 I = 1, N
219 D( I ) = DLAMDA( INDX( I ) )
220 Z( I ) = W( INDX( I ) )
221 40 CONTINUE
222 *
223 * Calculate the allowable deflation tolerance
224 *
225 IMAX = IDAMAX( N, Z, 1 )
226 JMAX = IDAMAX( N, D, 1 )
227 EPS = DLAMCH( 'Epsilon' )
228 TOL = EIGHT*EPS*ABS( D( JMAX ) )
229 *
230 * If the rank-1 modifier is small enough, no more needs to be done
231 * -- except to reorganize Q so that its columns correspond with the
232 * elements in D.
233 *
234 IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
235 K = 0
236 DO 50 J = 1, N
237 PERM( J ) = INDXQ( INDX( J ) )
238 CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
239 50 CONTINUE
240 CALL ZLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), LDQ )
241 RETURN
242 END IF
243 *
244 * If there are multiple eigenvalues then the problem deflates. Here
245 * the number of equal eigenvalues are found. As each equal
246 * eigenvalue is found, an elementary reflector is computed to rotate
247 * the corresponding eigensubspace so that the corresponding
248 * components of Z are zero in this new basis.
249 *
250 K = 0
251 K2 = N + 1
252 DO 60 J = 1, N
253 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
254 *
255 * Deflate due to small z component.
256 *
257 K2 = K2 - 1
258 INDXP( K2 ) = J
259 IF( J.EQ.N )
260 $ GO TO 100
261 ELSE
262 JLAM = J
263 GO TO 70
264 END IF
265 60 CONTINUE
266 70 CONTINUE
267 J = J + 1
268 IF( J.GT.N )
269 $ GO TO 90
270 IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
271 *
272 * Deflate due to small z component.
273 *
274 K2 = K2 - 1
275 INDXP( K2 ) = J
276 ELSE
277 *
278 * Check if eigenvalues are close enough to allow deflation.
279 *
280 S = Z( JLAM )
281 C = Z( J )
282 *
283 * Find sqrt(a**2+b**2) without overflow or
284 * destructive underflow.
285 *
286 TAU = DLAPY2( C, S )
287 T = D( J ) - D( JLAM )
288 C = C / TAU
289 S = -S / TAU
290 IF( ABS( T*C*S ).LE.TOL ) THEN
291 *
292 * Deflation is possible.
293 *
294 Z( J ) = TAU
295 Z( JLAM ) = ZERO
296 *
297 * Record the appropriate Givens rotation
298 *
299 GIVPTR = GIVPTR + 1
300 GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
301 GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
302 GIVNUM( 1, GIVPTR ) = C
303 GIVNUM( 2, GIVPTR ) = S
304 CALL ZDROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
305 $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
306 T = D( JLAM )*C*C + D( J )*S*S
307 D( J ) = D( JLAM )*S*S + D( J )*C*C
308 D( JLAM ) = T
309 K2 = K2 - 1
310 I = 1
311 80 CONTINUE
312 IF( K2+I.LE.N ) THEN
313 IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
314 INDXP( K2+I-1 ) = INDXP( K2+I )
315 INDXP( K2+I ) = JLAM
316 I = I + 1
317 GO TO 80
318 ELSE
319 INDXP( K2+I-1 ) = JLAM
320 END IF
321 ELSE
322 INDXP( K2+I-1 ) = JLAM
323 END IF
324 JLAM = J
325 ELSE
326 K = K + 1
327 W( K ) = Z( JLAM )
328 DLAMDA( K ) = D( JLAM )
329 INDXP( K ) = JLAM
330 JLAM = J
331 END IF
332 END IF
333 GO TO 70
334 90 CONTINUE
335 *
336 * Record the last eigenvalue.
337 *
338 K = K + 1
339 W( K ) = Z( JLAM )
340 DLAMDA( K ) = D( JLAM )
341 INDXP( K ) = JLAM
342 *
343 100 CONTINUE
344 *
345 * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
346 * and Q2 respectively. The eigenvalues/vectors which were not
347 * deflated go into the first K slots of DLAMDA and Q2 respectively,
348 * while those which were deflated go into the last N - K slots.
349 *
350 DO 110 J = 1, N
351 JP = INDXP( J )
352 DLAMDA( J ) = D( JP )
353 PERM( J ) = INDXQ( INDX( JP ) )
354 CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
355 110 CONTINUE
356 *
357 * The deflated eigenvalues and their corresponding vectors go back
358 * into the last N - K slots of D and Q respectively.
359 *
360 IF( K.LT.N ) THEN
361 CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
362 CALL ZLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ),
363 $ LDQ )
364 END IF
365 *
366 RETURN
367 *
368 * End of ZLAED8
369 *
370 END