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