1: SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W,
2: $ Q2, INDX, INDXC, INDXP, COLTYP, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, K, LDQ, N, N1
11: DOUBLE PRECISION RHO
12: * ..
13: * .. Array Arguments ..
14: INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ),
15: $ INDXQ( * )
16: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
17: $ W( * ), Z( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * DLAED2 merges the two sets of eigenvalues together into a single
24: * sorted set. Then it tries to deflate the size of the problem.
25: * There are two ways in which deflation can occur: when two or more
26: * eigenvalues are close together or if there is a tiny entry in the
27: * Z vector. For each such occurrence the order of the related secular
28: * equation problem is reduced by one.
29: *
30: * Arguments
31: * =========
32: *
33: * K (output) INTEGER
34: * The number of non-deflated eigenvalues, and the order of the
35: * related secular equation. 0 <= K <=N.
36: *
37: * N (input) INTEGER
38: * The dimension of the symmetric tridiagonal matrix. N >= 0.
39: *
40: * N1 (input) INTEGER
41: * The location of the last eigenvalue in the leading sub-matrix.
42: * min(1,N) <= N1 <= N/2.
43: *
44: * D (input/output) DOUBLE PRECISION array, dimension (N)
45: * On entry, D contains the eigenvalues of the two submatrices to
46: * be combined.
47: * On exit, D contains the trailing (N-K) updated eigenvalues
48: * (those which were deflated) sorted into increasing order.
49: *
50: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
51: * On entry, Q contains the eigenvectors of two submatrices in
52: * the two square blocks with corners at (1,1), (N1,N1)
53: * and (N1+1, N1+1), (N,N).
54: * On exit, Q contains the trailing (N-K) updated eigenvectors
55: * (those which were deflated) in its last N-K columns.
56: *
57: * LDQ (input) INTEGER
58: * The leading dimension of the array Q. LDQ >= max(1,N).
59: *
60: * INDXQ (input/output) INTEGER array, dimension (N)
61: * The permutation which separately sorts the two sub-problems
62: * in D into ascending order. Note that elements in the second
63: * half of this permutation must first have N1 added to their
64: * values. Destroyed on exit.
65: *
66: * RHO (input/output) DOUBLE PRECISION
67: * On entry, the off-diagonal element associated with the rank-1
68: * cut which originally split the two submatrices which are now
69: * being recombined.
70: * On exit, RHO has been modified to the value required by
71: * DLAED3.
72: *
73: * Z (input) DOUBLE PRECISION array, dimension (N)
74: * On entry, Z 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).
77: * On exit, the contents of Z have been destroyed by the updating
78: * process.
79: *
80: * DLAMDA (output) DOUBLE PRECISION array, dimension (N)
81: * A copy of the first K eigenvalues which will be used by
82: * DLAED3 to form the secular equation.
83: *
84: * W (output) DOUBLE PRECISION array, dimension (N)
85: * The first k values of the final deflation-altered z-vector
86: * which will be passed to DLAED3.
87: *
88: * Q2 (output) DOUBLE PRECISION array, dimension (N1**2+(N-N1)**2)
89: * A copy of the first K eigenvectors which will be used by
90: * DLAED3 in a matrix multiply (DGEMM) to solve for the new
91: * eigenvectors.
92: *
93: * INDX (workspace) INTEGER array, dimension (N)
94: * The permutation used to sort the contents of DLAMDA into
95: * ascending order.
96: *
97: * INDXC (output) INTEGER array, dimension (N)
98: * The permutation used to arrange the columns of the deflated
99: * Q matrix into three groups: the first group contains non-zero
100: * elements only at and above N1, the second contains
101: * non-zero elements only below N1, and the third is dense.
102: *
103: * INDXP (workspace) INTEGER array, dimension (N)
104: * The permutation used to place deflated values of D at the end
105: * of the array. INDXP(1:K) points to the nondeflated D-values
106: * and INDXP(K+1:N) points to the deflated eigenvalues.
107: *
108: * COLTYP (workspace/output) INTEGER array, dimension (N)
109: * During execution, a label which will indicate which of the
110: * following types a column in the Q2 matrix is:
111: * 1 : non-zero in the upper half only;
112: * 2 : dense;
113: * 3 : non-zero in the lower half only;
114: * 4 : deflated.
115: * On exit, COLTYP(i) is the number of columns of type i,
116: * for i=1 to 4 only.
117: *
118: * INFO (output) INTEGER
119: * = 0: successful exit.
120: * < 0: if INFO = -i, the i-th argument had an illegal value.
121: *
122: * Further Details
123: * ===============
124: *
125: * Based on contributions by
126: * Jeff Rutter, Computer Science Division, University of California
127: * at Berkeley, USA
128: * Modified by Francoise Tisseur, University of Tennessee.
129: *
130: * =====================================================================
131: *
132: * .. Parameters ..
133: DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
134: PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
135: $ TWO = 2.0D0, EIGHT = 8.0D0 )
136: * ..
137: * .. Local Arrays ..
138: INTEGER CTOT( 4 ), PSM( 4 )
139: * ..
140: * .. Local Scalars ..
141: INTEGER CT, I, IMAX, IQ1, IQ2, J, JMAX, JS, K2, N1P1,
142: $ N2, NJ, PJ
143: DOUBLE PRECISION C, EPS, S, T, TAU, TOL
144: * ..
145: * .. External Functions ..
146: INTEGER IDAMAX
147: DOUBLE PRECISION DLAMCH, DLAPY2
148: EXTERNAL IDAMAX, DLAMCH, DLAPY2
149: * ..
150: * .. External Subroutines ..
151: EXTERNAL DCOPY, DLACPY, DLAMRG, DROT, DSCAL, XERBLA
152: * ..
153: * .. Intrinsic Functions ..
154: INTRINSIC ABS, MAX, MIN, SQRT
155: * ..
156: * .. Executable Statements ..
157: *
158: * Test the input parameters.
159: *
160: INFO = 0
161: *
162: IF( N.LT.0 ) THEN
163: INFO = -2
164: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
165: INFO = -6
166: ELSE IF( MIN( 1, ( N / 2 ) ).GT.N1 .OR. ( N / 2 ).LT.N1 ) THEN
167: INFO = -3
168: END IF
169: IF( INFO.NE.0 ) THEN
170: CALL XERBLA( 'DLAED2', -INFO )
171: RETURN
172: END IF
173: *
174: * Quick return if possible
175: *
176: IF( N.EQ.0 )
177: $ RETURN
178: *
179: N2 = N - N1
180: N1P1 = N1 + 1
181: *
182: IF( RHO.LT.ZERO ) THEN
183: CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
184: END IF
185: *
186: * Normalize z so that norm(z) = 1. Since z is the concatenation of
187: * two normalized vectors, norm2(z) = sqrt(2).
188: *
189: T = ONE / SQRT( TWO )
190: CALL DSCAL( N, T, Z, 1 )
191: *
192: * RHO = ABS( norm(z)**2 * RHO )
193: *
194: RHO = ABS( TWO*RHO )
195: *
196: * Sort the eigenvalues into increasing order
197: *
198: DO 10 I = N1P1, N
199: INDXQ( I ) = INDXQ( I ) + N1
200: 10 CONTINUE
201: *
202: * re-integrate the deflated parts from the last pass
203: *
204: DO 20 I = 1, N
205: DLAMDA( I ) = D( INDXQ( I ) )
206: 20 CONTINUE
207: CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDXC )
208: DO 30 I = 1, N
209: INDX( I ) = INDXQ( INDXC( I ) )
210: 30 CONTINUE
211: *
212: * Calculate the allowable deflation tolerance
213: *
214: IMAX = IDAMAX( N, Z, 1 )
215: JMAX = IDAMAX( N, D, 1 )
216: EPS = DLAMCH( 'Epsilon' )
217: TOL = EIGHT*EPS*MAX( ABS( D( JMAX ) ), ABS( Z( IMAX ) ) )
218: *
219: * If the rank-1 modifier is small enough, no more needs to be done
220: * except to reorganize Q so that its columns correspond with the
221: * elements in D.
222: *
223: IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
224: K = 0
225: IQ2 = 1
226: DO 40 J = 1, N
227: I = INDX( J )
228: CALL DCOPY( N, Q( 1, I ), 1, Q2( IQ2 ), 1 )
229: DLAMDA( J ) = D( I )
230: IQ2 = IQ2 + N
231: 40 CONTINUE
232: CALL DLACPY( 'A', N, N, Q2, N, Q, LDQ )
233: CALL DCOPY( N, DLAMDA, 1, D, 1 )
234: GO TO 190
235: END IF
236: *
237: * If there are multiple eigenvalues then the problem deflates. Here
238: * the number of equal eigenvalues are found. As each equal
239: * eigenvalue is found, an elementary reflector is computed to rotate
240: * the corresponding eigensubspace so that the corresponding
241: * components of Z are zero in this new basis.
242: *
243: DO 50 I = 1, N1
244: COLTYP( I ) = 1
245: 50 CONTINUE
246: DO 60 I = N1P1, N
247: COLTYP( I ) = 3
248: 60 CONTINUE
249: *
250: *
251: K = 0
252: K2 = N + 1
253: DO 70 J = 1, N
254: NJ = INDX( J )
255: IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
256: *
257: * Deflate due to small z component.
258: *
259: K2 = K2 - 1
260: COLTYP( NJ ) = 4
261: INDXP( K2 ) = NJ
262: IF( J.EQ.N )
263: $ GO TO 100
264: ELSE
265: PJ = NJ
266: GO TO 80
267: END IF
268: 70 CONTINUE
269: 80 CONTINUE
270: J = J + 1
271: NJ = INDX( J )
272: IF( J.GT.N )
273: $ GO TO 100
274: IF( RHO*ABS( Z( NJ ) ).LE.TOL ) THEN
275: *
276: * Deflate due to small z component.
277: *
278: K2 = K2 - 1
279: COLTYP( NJ ) = 4
280: INDXP( K2 ) = NJ
281: ELSE
282: *
283: * Check if eigenvalues are close enough to allow deflation.
284: *
285: S = Z( PJ )
286: C = Z( NJ )
287: *
288: * Find sqrt(a**2+b**2) without overflow or
289: * destructive underflow.
290: *
291: TAU = DLAPY2( C, S )
292: T = D( NJ ) - D( PJ )
293: C = C / TAU
294: S = -S / TAU
295: IF( ABS( T*C*S ).LE.TOL ) THEN
296: *
297: * Deflation is possible.
298: *
299: Z( NJ ) = TAU
300: Z( PJ ) = ZERO
301: IF( COLTYP( NJ ).NE.COLTYP( PJ ) )
302: $ COLTYP( NJ ) = 2
303: COLTYP( PJ ) = 4
304: CALL DROT( N, Q( 1, PJ ), 1, Q( 1, NJ ), 1, C, S )
305: T = D( PJ )*C**2 + D( NJ )*S**2
306: D( NJ ) = D( PJ )*S**2 + D( NJ )*C**2
307: D( PJ ) = T
308: K2 = K2 - 1
309: I = 1
310: 90 CONTINUE
311: IF( K2+I.LE.N ) THEN
312: IF( D( PJ ).LT.D( INDXP( K2+I ) ) ) THEN
313: INDXP( K2+I-1 ) = INDXP( K2+I )
314: INDXP( K2+I ) = PJ
315: I = I + 1
316: GO TO 90
317: ELSE
318: INDXP( K2+I-1 ) = PJ
319: END IF
320: ELSE
321: INDXP( K2+I-1 ) = PJ
322: END IF
323: PJ = NJ
324: ELSE
325: K = K + 1
326: DLAMDA( K ) = D( PJ )
327: W( K ) = Z( PJ )
328: INDXP( K ) = PJ
329: PJ = NJ
330: END IF
331: END IF
332: GO TO 80
333: 100 CONTINUE
334: *
335: * Record the last eigenvalue.
336: *
337: K = K + 1
338: DLAMDA( K ) = D( PJ )
339: W( K ) = Z( PJ )
340: INDXP( K ) = PJ
341: *
342: * Count up the total number of the various types of columns, then
343: * form a permutation which positions the four column types into
344: * four uniform groups (although one or more of these groups may be
345: * empty).
346: *
347: DO 110 J = 1, 4
348: CTOT( J ) = 0
349: 110 CONTINUE
350: DO 120 J = 1, N
351: CT = COLTYP( J )
352: CTOT( CT ) = CTOT( CT ) + 1
353: 120 CONTINUE
354: *
355: * PSM(*) = Position in SubMatrix (of types 1 through 4)
356: *
357: PSM( 1 ) = 1
358: PSM( 2 ) = 1 + CTOT( 1 )
359: PSM( 3 ) = PSM( 2 ) + CTOT( 2 )
360: PSM( 4 ) = PSM( 3 ) + CTOT( 3 )
361: K = N - CTOT( 4 )
362: *
363: * Fill out the INDXC array so that the permutation which it induces
364: * will place all type-1 columns first, all type-2 columns next,
365: * then all type-3's, and finally all type-4's.
366: *
367: DO 130 J = 1, N
368: JS = INDXP( J )
369: CT = COLTYP( JS )
370: INDX( PSM( CT ) ) = JS
371: INDXC( PSM( CT ) ) = J
372: PSM( CT ) = PSM( CT ) + 1
373: 130 CONTINUE
374: *
375: * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
376: * and Q2 respectively. The eigenvalues/vectors which were not
377: * deflated go into the first K slots of DLAMDA and Q2 respectively,
378: * while those which were deflated go into the last N - K slots.
379: *
380: I = 1
381: IQ1 = 1
382: IQ2 = 1 + ( CTOT( 1 )+CTOT( 2 ) )*N1
383: DO 140 J = 1, CTOT( 1 )
384: JS = INDX( I )
385: CALL DCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 )
386: Z( I ) = D( JS )
387: I = I + 1
388: IQ1 = IQ1 + N1
389: 140 CONTINUE
390: *
391: DO 150 J = 1, CTOT( 2 )
392: JS = INDX( I )
393: CALL DCOPY( N1, Q( 1, JS ), 1, Q2( IQ1 ), 1 )
394: CALL DCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 )
395: Z( I ) = D( JS )
396: I = I + 1
397: IQ1 = IQ1 + N1
398: IQ2 = IQ2 + N2
399: 150 CONTINUE
400: *
401: DO 160 J = 1, CTOT( 3 )
402: JS = INDX( I )
403: CALL DCOPY( N2, Q( N1+1, JS ), 1, Q2( IQ2 ), 1 )
404: Z( I ) = D( JS )
405: I = I + 1
406: IQ2 = IQ2 + N2
407: 160 CONTINUE
408: *
409: IQ1 = IQ2
410: DO 170 J = 1, CTOT( 4 )
411: JS = INDX( I )
412: CALL DCOPY( N, Q( 1, JS ), 1, Q2( IQ2 ), 1 )
413: IQ2 = IQ2 + N
414: Z( I ) = D( JS )
415: I = I + 1
416: 170 CONTINUE
417: *
418: * The deflated eigenvalues and their corresponding vectors go back
419: * into the last N - K slots of D and Q respectively.
420: *
421: CALL DLACPY( 'A', N, CTOT( 4 ), Q2( IQ1 ), N, Q( 1, K+1 ), LDQ )
422: CALL DCOPY( N-K, Z( K+1 ), 1, D( K+1 ), 1 )
423: *
424: * Copy CTOT into COLTYP for referencing in DLAED3.
425: *
426: DO 180 J = 1, 4
427: COLTYP( J ) = CTOT( J )
428: 180 CONTINUE
429: *
430: 190 CONTINUE
431: RETURN
432: *
433: * End of DLAED2
434: *
435: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack