1: *> \brief \b ZLAED8 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 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: *> \date December 2016
223: *
224: *> \ingroup complex16OTHERcomputational
225: *
226: * =====================================================================
227: SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA,
228: $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR,
229: $ GIVCOL, GIVNUM, INFO )
230: *
231: * -- LAPACK computational routine (version 3.7.0) --
232: * -- LAPACK is a software package provided by Univ. of Tennessee, --
233: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234: * December 2016
235: *
236: * .. Scalar Arguments ..
237: INTEGER CUTPNT, GIVPTR, INFO, K, LDQ, LDQ2, N, QSIZ
238: DOUBLE PRECISION RHO
239: * ..
240: * .. Array Arguments ..
241: INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ),
242: $ INDXQ( * ), PERM( * )
243: DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ),
244: $ Z( * )
245: COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * )
246: * ..
247: *
248: * =====================================================================
249: *
250: * .. Parameters ..
251: DOUBLE PRECISION MONE, ZERO, ONE, TWO, EIGHT
252: PARAMETER ( MONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0,
253: $ TWO = 2.0D0, EIGHT = 8.0D0 )
254: * ..
255: * .. Local Scalars ..
256: INTEGER I, IMAX, J, JLAM, JMAX, JP, K2, N1, N1P1, N2
257: DOUBLE PRECISION C, EPS, S, T, TAU, TOL
258: * ..
259: * .. External Functions ..
260: INTEGER IDAMAX
261: DOUBLE PRECISION DLAMCH, DLAPY2
262: EXTERNAL IDAMAX, DLAMCH, DLAPY2
263: * ..
264: * .. External Subroutines ..
265: EXTERNAL DCOPY, DLAMRG, DSCAL, XERBLA, ZCOPY, ZDROT,
266: $ ZLACPY
267: * ..
268: * .. Intrinsic Functions ..
269: INTRINSIC ABS, MAX, MIN, SQRT
270: * ..
271: * .. Executable Statements ..
272: *
273: * Test the input parameters.
274: *
275: INFO = 0
276: *
277: IF( N.LT.0 ) THEN
278: INFO = -2
279: ELSE IF( QSIZ.LT.N ) THEN
280: INFO = -3
281: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
282: INFO = -5
283: ELSE IF( CUTPNT.LT.MIN( 1, N ) .OR. CUTPNT.GT.N ) THEN
284: INFO = -8
285: ELSE IF( LDQ2.LT.MAX( 1, N ) ) THEN
286: INFO = -12
287: END IF
288: IF( INFO.NE.0 ) THEN
289: CALL XERBLA( 'ZLAED8', -INFO )
290: RETURN
291: END IF
292: *
293: * Need to initialize GIVPTR to O here in case of quick exit
294: * to prevent an unspecified code behavior (usually sigfault)
295: * when IWORK array on entry to *stedc is not zeroed
296: * (or at least some IWORK entries which used in *laed7 for GIVPTR).
297: *
298: GIVPTR = 0
299: *
300: * Quick return if possible
301: *
302: IF( N.EQ.0 )
303: $ RETURN
304: *
305: N1 = CUTPNT
306: N2 = N - N1
307: N1P1 = N1 + 1
308: *
309: IF( RHO.LT.ZERO ) THEN
310: CALL DSCAL( N2, MONE, Z( N1P1 ), 1 )
311: END IF
312: *
313: * Normalize z so that norm(z) = 1
314: *
315: T = ONE / SQRT( TWO )
316: DO 10 J = 1, N
317: INDX( J ) = J
318: 10 CONTINUE
319: CALL DSCAL( N, T, Z, 1 )
320: RHO = ABS( TWO*RHO )
321: *
322: * Sort the eigenvalues into increasing order
323: *
324: DO 20 I = CUTPNT + 1, N
325: INDXQ( I ) = INDXQ( I ) + CUTPNT
326: 20 CONTINUE
327: DO 30 I = 1, N
328: DLAMDA( I ) = D( INDXQ( I ) )
329: W( I ) = Z( INDXQ( I ) )
330: 30 CONTINUE
331: I = 1
332: J = CUTPNT + 1
333: CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX )
334: DO 40 I = 1, N
335: D( I ) = DLAMDA( INDX( I ) )
336: Z( I ) = W( INDX( I ) )
337: 40 CONTINUE
338: *
339: * Calculate the allowable deflation tolerance
340: *
341: IMAX = IDAMAX( N, Z, 1 )
342: JMAX = IDAMAX( N, D, 1 )
343: EPS = DLAMCH( 'Epsilon' )
344: TOL = EIGHT*EPS*ABS( D( JMAX ) )
345: *
346: * If the rank-1 modifier is small enough, no more needs to be done
347: * -- except to reorganize Q so that its columns correspond with the
348: * elements in D.
349: *
350: IF( RHO*ABS( Z( IMAX ) ).LE.TOL ) THEN
351: K = 0
352: DO 50 J = 1, N
353: PERM( J ) = INDXQ( INDX( J ) )
354: CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
355: 50 CONTINUE
356: CALL ZLACPY( 'A', QSIZ, N, Q2( 1, 1 ), LDQ2, Q( 1, 1 ), LDQ )
357: RETURN
358: END IF
359: *
360: * If there are multiple eigenvalues then the problem deflates. Here
361: * the number of equal eigenvalues are found. As each equal
362: * eigenvalue is found, an elementary reflector is computed to rotate
363: * the corresponding eigensubspace so that the corresponding
364: * components of Z are zero in this new basis.
365: *
366: K = 0
367: K2 = N + 1
368: DO 60 J = 1, N
369: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
370: *
371: * Deflate due to small z component.
372: *
373: K2 = K2 - 1
374: INDXP( K2 ) = J
375: IF( J.EQ.N )
376: $ GO TO 100
377: ELSE
378: JLAM = J
379: GO TO 70
380: END IF
381: 60 CONTINUE
382: 70 CONTINUE
383: J = J + 1
384: IF( J.GT.N )
385: $ GO TO 90
386: IF( RHO*ABS( Z( J ) ).LE.TOL ) THEN
387: *
388: * Deflate due to small z component.
389: *
390: K2 = K2 - 1
391: INDXP( K2 ) = J
392: ELSE
393: *
394: * Check if eigenvalues are close enough to allow deflation.
395: *
396: S = Z( JLAM )
397: C = Z( J )
398: *
399: * Find sqrt(a**2+b**2) without overflow or
400: * destructive underflow.
401: *
402: TAU = DLAPY2( C, S )
403: T = D( J ) - D( JLAM )
404: C = C / TAU
405: S = -S / TAU
406: IF( ABS( T*C*S ).LE.TOL ) THEN
407: *
408: * Deflation is possible.
409: *
410: Z( J ) = TAU
411: Z( JLAM ) = ZERO
412: *
413: * Record the appropriate Givens rotation
414: *
415: GIVPTR = GIVPTR + 1
416: GIVCOL( 1, GIVPTR ) = INDXQ( INDX( JLAM ) )
417: GIVCOL( 2, GIVPTR ) = INDXQ( INDX( J ) )
418: GIVNUM( 1, GIVPTR ) = C
419: GIVNUM( 2, GIVPTR ) = S
420: CALL ZDROT( QSIZ, Q( 1, INDXQ( INDX( JLAM ) ) ), 1,
421: $ Q( 1, INDXQ( INDX( J ) ) ), 1, C, S )
422: T = D( JLAM )*C*C + D( J )*S*S
423: D( J ) = D( JLAM )*S*S + D( J )*C*C
424: D( JLAM ) = T
425: K2 = K2 - 1
426: I = 1
427: 80 CONTINUE
428: IF( K2+I.LE.N ) THEN
429: IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
430: INDXP( K2+I-1 ) = INDXP( K2+I )
431: INDXP( K2+I ) = JLAM
432: I = I + 1
433: GO TO 80
434: ELSE
435: INDXP( K2+I-1 ) = JLAM
436: END IF
437: ELSE
438: INDXP( K2+I-1 ) = JLAM
439: END IF
440: JLAM = J
441: ELSE
442: K = K + 1
443: W( K ) = Z( JLAM )
444: DLAMDA( K ) = D( JLAM )
445: INDXP( K ) = JLAM
446: JLAM = J
447: END IF
448: END IF
449: GO TO 70
450: 90 CONTINUE
451: *
452: * Record the last eigenvalue.
453: *
454: K = K + 1
455: W( K ) = Z( JLAM )
456: DLAMDA( K ) = D( JLAM )
457: INDXP( K ) = JLAM
458: *
459: 100 CONTINUE
460: *
461: * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
462: * and Q2 respectively. The eigenvalues/vectors which were not
463: * deflated go into the first K slots of DLAMDA and Q2 respectively,
464: * while those which were deflated go into the last N - K slots.
465: *
466: DO 110 J = 1, N
467: JP = INDXP( J )
468: DLAMDA( J ) = D( JP )
469: PERM( J ) = INDXQ( INDX( JP ) )
470: CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 )
471: 110 CONTINUE
472: *
473: * The deflated eigenvalues and their corresponding vectors go back
474: * into the last N - K slots of D and Q respectively.
475: *
476: IF( K.LT.N ) THEN
477: CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 )
478: CALL ZLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ),
479: $ LDQ )
480: END IF
481: *
482: RETURN
483: *
484: * End of ZLAED8
485: *
486: 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