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