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