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.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, 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: * Need to initialize GIVPTR to O here in case of quick exit
193: * to prevent an unspecified code behavior (usually sigfault)
194: * when IWORK array on entry to *stedc is not zeroed
195: * (or at least some IWORK entries which used in *laed7 for GIVPTR).
196: *
197: GIVPTR = 0
198: *
199: * Quick return if possible
200: *
201: IF( N.EQ.0 )
202: $ RETURN
203: *
204: N1 = CUTPNT
205: N2 = N - N1
206: N1P1 = N1 + 1
207: *
208: IF( RHO.LT.ZERO ) THEN
209: CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
210: END IF
211: *
212: * Normalize z so that norm(z) = 1
213: *
214: T = ONE / SQRT( TWO )
215: DO 10 J = 1, N
216: INDX( J ) = J
217: 10 CONTINUE
218: CALL DSCAL( N, T, Z, 1 )
219: RHO = ABS( TWO*RHO )
220: *
221: * Sort the eigenvalues into increasing order
222: *
223: DO 20 I = CUTPNT + 1, N
224: INDXQ( I ) = INDXQ( I ) + CUTPNT
225: 20 CONTINUE
226: DO 30 I = 1, N
227: DLAMDA( I ) = D( INDXQ( I ) )
228: W( I ) = Z( INDXQ( I ) )
229: 30 CONTINUE
230: I = 1
231: J = CUTPNT + 1
232: CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
233: DO 40 I = 1, N
234: D( I ) = DLAMDA( INDX( I ) )
235: Z( I ) = W( INDX( I ) )
236: 40 CONTINUE
237: *
238: * Calculate the allowable deflation tolerence
239: *
240: IMAX = IDAMAX( N, Z, 1 )
241: JMAX = IDAMAX( N, D, 1 )
242: EPS = DLAMCH( 'Epsilon' )
243: TOL = EIGHT*EPS*ABS( D( JMAX ) )
244: *
245: * If the rank-1 modifier is small enough, no more needs to be done
246: * except to reorganize Q so that its columns correspond with the
247: * elements in D.
248: *
249: IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
250: K = 0
251: IF( ICOMPQ.EQ.0 ) THEN
252: DO 50 J = 1, N
253: PERM( J ) = INDXQ( INDX( J ) )
254: 50 CONTINUE
255: ELSE
256: DO 60 J = 1, N
257: PERM( J ) = INDXQ( INDX( J ) )
258: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
259: 60 CONTINUE
260: CALL DLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ),
261: $ LDQ )
262: END IF
263: RETURN
264: END IF
265: *
266: * If there are multiple eigenvalues then the problem deflates. Here
267: * the number of equal eigenvalues are found. As each equal
268: * eigenvalue is found, an elementary reflector is computed to rotate
269: * the corresponding eigensubspace so that the corresponding
270: * components of Z are zero in this new basis.
271: *
272: K = 0
273: K2 = N + 1
274: DO 70 J = 1, N
275: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
276: *
277: * Deflate due to small z component.
278: *
279: K2 = K2 - 1
280: INDXP( K2 ) = J
281: IF( J.EQ.N )
282: $ GO TO 110
283: ELSE
284: JLAM = J
285: GO TO 80
286: END IF
287: 70 CONTINUE
288: 80 CONTINUE
289: J = J + 1
290: IF( J.GT.N )
291: $ GO TO 100
292: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
293: *
294: * Deflate due to small z component.
295: *
296: K2 = K2 - 1
297: INDXP( K2 ) = J
298: ELSE
299: *
300: * Check if eigenvalues are close enough to allow deflation.
301: *
302: S = Z( JLAM )
303: C = Z( J )
304: *
305: * Find sqrt(a**2+b**2) without overflow or
306: * destructive underflow.
307: *
308: TAU = DLAPY2( C, S )
309: T = D( J ) - D( JLAM )
310: C = C / TAU
311: S = -S / TAU
312: IF( ABS( T*C*S ).LE.TOL ) THEN
313: *
314: * Deflation is possible.
315: *
316: Z( J ) = TAU
317: Z( JLAM ) = ZERO
318: *
319: * Record the appropriate Givens rotation
320: *
321: GIVPTR = GIVPTR + 1
322: GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
323: GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
324: GIVNUM( 1, GIVPTR ) = C
325: GIVNUM( 2, GIVPTR ) = S
326: IF( ICOMPQ.EQ.1 ) THEN
327: CALL DROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
328: $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
329: END IF
330: T = D( JLAM )*C*C + D( J )*S*S
331: D( J ) = D( JLAM )*S*S + D( J )*C*C
332: D( JLAM ) = T
333: K2 = K2 - 1
334: I = 1
335: 90 CONTINUE
336: IF( K2+I.LE.N ) THEN
337: IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
338: INDXP( K2+I-1 ) = INDXP( K2+I )
339: INDXP( K2+I ) = JLAM
340: I = I + 1
341: GO TO 90
342: ELSE
343: INDXP( K2+I-1 ) = JLAM
344: END IF
345: ELSE
346: INDXP( K2+I-1 ) = JLAM
347: END IF
348: JLAM = J
349: ELSE
350: K = K + 1
351: W( K ) = Z( JLAM )
352: DLAMDA( K ) = D( JLAM )
353: INDXP( K ) = JLAM
354: JLAM = J
355: END IF
356: END IF
357: GO TO 80
358: 100 CONTINUE
359: *
360: * Record the last eigenvalue.
361: *
362: K = K + 1
363: W( K ) = Z( JLAM )
364: DLAMDA( K ) = D( JLAM )
365: INDXP( K ) = JLAM
366: *
367: 110 CONTINUE
368: *
369: * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
370: * and Q2 respectively. The eigenvalues/vectors which were not
371: * deflated go into the first K slots of DLAMDA and Q2 respectively,
372: * while those which were deflated go into the last N - K slots.
373: *
374: IF( ICOMPQ.EQ.0 ) THEN
375: DO 120 J = 1, N
376: JP = INDXP( J )
377: DLAMDA( J ) = D( JP )
378: PERM( J ) = INDXQ( INDX( JP ) )
379: 120 CONTINUE
380: ELSE
381: DO 130 J = 1, N
382: JP = INDXP( J )
383: DLAMDA( J ) = D( JP )
384: PERM( J ) = INDXQ( INDX( JP ) )
385: CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
386: 130 CONTINUE
387: END IF
388: *
389: * The deflated eigenvalues and their corresponding vectors go back
390: * into the last N - K slots of D and Q respectively.
391: *
392: IF( K.LT.N ) THEN
393: IF( ICOMPQ.EQ.0 ) THEN
394: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
395: ELSE
396: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
397: CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2,
398: $ Q( 1, K+1 ), LDQ )
399: END IF
400: END IF
401: *
402: RETURN
403: *
404: * End of DLAED8
405: *
406: 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