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