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