Annotation of rpl/lapack/lapack/zlaed7.f, revision 1.20
1.20 ! bertrand 1: *> \brief \b ZLAED7 used by ZSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a rank-one symmetric matrix. Used when the original matrix is dense.
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 9: *> Download ZLAED7 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed7.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed7.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed7.f">
1.9 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
22: * LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
23: * GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
24: * INFO )
1.17 bertrand 25: *
1.9 bertrand 26: * .. Scalar Arguments ..
27: * INTEGER CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
28: * $ TLVLS
29: * DOUBLE PRECISION RHO
30: * ..
31: * .. Array Arguments ..
32: * INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
33: * $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
34: * DOUBLE PRECISION D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
35: * COMPLEX*16 Q( LDQ, * ), WORK( * )
36: * ..
1.17 bertrand 37: *
1.9 bertrand 38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> ZLAED7 computes the updated eigensystem of a diagonal
45: *> matrix after modification by a rank-one symmetric matrix. This
46: *> routine is used only for the eigenproblem which requires all
47: *> eigenvalues and optionally eigenvectors of a dense or banded
48: *> Hermitian matrix that has been reduced to tridiagonal form.
49: *>
50: *> T = Q(in) ( D(in) + RHO * Z*Z**H ) Q**H(in) = Q(out) * D(out) * Q**H(out)
51: *>
52: *> where Z = Q**Hu, u is a vector of length N with ones in the
53: *> CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
54: *>
55: *> The eigenvectors of the original matrix are stored in Q, and the
56: *> eigenvalues are in D. The algorithm consists of three stages:
57: *>
58: *> The first stage consists of deflating the size of the problem
59: *> when there are multiple eigenvalues or if there is a zero in
1.15 bertrand 60: *> the Z vector. For each such occurrence the dimension of the
1.9 bertrand 61: *> secular equation problem is reduced by one. This stage is
62: *> performed by the routine DLAED2.
63: *>
64: *> The second stage consists of calculating the updated
65: *> eigenvalues. This is done by finding the roots of the secular
66: *> equation via the routine DLAED4 (as called by SLAED3).
67: *> This routine also calculates the eigenvectors of the current
68: *> problem.
69: *>
70: *> The final stage consists of computing the updated eigenvectors
71: *> directly using the updated eigenvalues. The eigenvectors for
72: *> the current problem are multiplied with the eigenvectors from
73: *> the overall problem.
74: *> \endverbatim
75: *
76: * Arguments:
77: * ==========
78: *
79: *> \param[in] N
80: *> \verbatim
81: *> N is INTEGER
82: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
83: *> \endverbatim
84: *>
85: *> \param[in] CUTPNT
86: *> \verbatim
87: *> CUTPNT is INTEGER
88: *> Contains the location of the last eigenvalue in the leading
89: *> sub-matrix. min(1,N) <= CUTPNT <= N.
90: *> \endverbatim
91: *>
92: *> \param[in] QSIZ
93: *> \verbatim
94: *> QSIZ is INTEGER
95: *> The dimension of the unitary matrix used to reduce
96: *> the full matrix to tridiagonal form. QSIZ >= N.
97: *> \endverbatim
98: *>
99: *> \param[in] TLVLS
100: *> \verbatim
101: *> TLVLS is INTEGER
102: *> The total number of merging levels in the overall divide and
103: *> conquer tree.
104: *> \endverbatim
105: *>
106: *> \param[in] CURLVL
107: *> \verbatim
108: *> CURLVL is INTEGER
109: *> The current level in the overall merge routine,
110: *> 0 <= curlvl <= tlvls.
111: *> \endverbatim
112: *>
113: *> \param[in] CURPBM
114: *> \verbatim
115: *> CURPBM is INTEGER
116: *> The current problem in the current level in the overall
117: *> merge routine (counting from upper left to lower right).
118: *> \endverbatim
119: *>
120: *> \param[in,out] D
121: *> \verbatim
122: *> D is DOUBLE PRECISION array, dimension (N)
123: *> On entry, the eigenvalues of the rank-1-perturbed matrix.
124: *> On exit, the eigenvalues of the repaired matrix.
125: *> \endverbatim
126: *>
127: *> \param[in,out] Q
128: *> \verbatim
129: *> Q is COMPLEX*16 array, dimension (LDQ,N)
130: *> On entry, the eigenvectors of the rank-1-perturbed matrix.
131: *> On exit, the eigenvectors of the repaired tridiagonal matrix.
132: *> \endverbatim
133: *>
134: *> \param[in] LDQ
135: *> \verbatim
136: *> LDQ is INTEGER
137: *> The leading dimension of the array Q. LDQ >= max(1,N).
138: *> \endverbatim
139: *>
140: *> \param[in] RHO
141: *> \verbatim
142: *> RHO is DOUBLE PRECISION
143: *> Contains the subdiagonal element used to create the rank-1
144: *> modification.
145: *> \endverbatim
146: *>
147: *> \param[out] INDXQ
148: *> \verbatim
149: *> INDXQ is INTEGER array, dimension (N)
150: *> This contains the permutation which will reintegrate the
151: *> subproblem just solved back into sorted order,
152: *> ie. D( INDXQ( I = 1, N ) ) will be in ascending order.
153: *> \endverbatim
154: *>
155: *> \param[out] IWORK
156: *> \verbatim
157: *> IWORK is INTEGER array, dimension (4*N)
158: *> \endverbatim
159: *>
160: *> \param[out] RWORK
161: *> \verbatim
162: *> RWORK is DOUBLE PRECISION array,
163: *> dimension (3*N+2*QSIZ*N)
164: *> \endverbatim
165: *>
166: *> \param[out] WORK
167: *> \verbatim
168: *> WORK is COMPLEX*16 array, dimension (QSIZ*N)
169: *> \endverbatim
170: *>
171: *> \param[in,out] QSTORE
172: *> \verbatim
173: *> QSTORE is DOUBLE PRECISION array, dimension (N**2+1)
174: *> Stores eigenvectors of submatrices encountered during
175: *> divide and conquer, packed together. QPTR points to
176: *> beginning of the submatrices.
177: *> \endverbatim
178: *>
179: *> \param[in,out] QPTR
180: *> \verbatim
181: *> QPTR is INTEGER array, dimension (N+2)
182: *> List of indices pointing to beginning of submatrices stored
183: *> in QSTORE. The submatrices are numbered starting at the
184: *> bottom left of the divide and conquer tree, from left to
185: *> right and bottom to top.
186: *> \endverbatim
187: *>
188: *> \param[in] PRMPTR
189: *> \verbatim
190: *> PRMPTR is INTEGER array, dimension (N lg N)
191: *> Contains a list of pointers which indicate where in PERM a
192: *> level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
193: *> indicates the size of the permutation and also the size of
194: *> the full, non-deflated problem.
195: *> \endverbatim
196: *>
197: *> \param[in] PERM
198: *> \verbatim
199: *> PERM is INTEGER array, dimension (N lg N)
200: *> Contains the permutations (from deflation and sorting) to be
201: *> applied to each eigenblock.
202: *> \endverbatim
203: *>
204: *> \param[in] GIVPTR
205: *> \verbatim
206: *> GIVPTR is INTEGER array, dimension (N lg N)
207: *> Contains a list of pointers which indicate where in GIVCOL a
208: *> level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
209: *> indicates the number of Givens rotations.
210: *> \endverbatim
211: *>
212: *> \param[in] GIVCOL
213: *> \verbatim
214: *> GIVCOL is INTEGER array, dimension (2, N lg N)
215: *> Each pair of numbers indicates a pair of columns to take place
216: *> in a Givens rotation.
217: *> \endverbatim
218: *>
219: *> \param[in] GIVNUM
220: *> \verbatim
221: *> GIVNUM is DOUBLE PRECISION array, dimension (2, N lg N)
222: *> Each number indicates the S value to be used in the
223: *> corresponding Givens rotation.
224: *> \endverbatim
225: *>
226: *> \param[out] INFO
227: *> \verbatim
228: *> INFO is INTEGER
229: *> = 0: successful exit.
230: *> < 0: if INFO = -i, the i-th argument had an illegal value.
231: *> > 0: if INFO = 1, an eigenvalue did not converge
232: *> \endverbatim
233: *
234: * Authors:
235: * ========
236: *
1.17 bertrand 237: *> \author Univ. of Tennessee
238: *> \author Univ. of California Berkeley
239: *> \author Univ. of Colorado Denver
240: *> \author NAG Ltd.
1.9 bertrand 241: *
242: *> \ingroup complex16OTHERcomputational
243: *
244: * =====================================================================
1.1 bertrand 245: SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
246: $ LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
247: $ GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
248: $ INFO )
249: *
1.20 ! bertrand 250: * -- LAPACK computational routine --
1.1 bertrand 251: * -- LAPACK is a software package provided by Univ. of Tennessee, --
252: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253: *
254: * .. Scalar Arguments ..
255: INTEGER CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
256: $ TLVLS
257: DOUBLE PRECISION RHO
258: * ..
259: * .. Array Arguments ..
260: INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
261: $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
262: DOUBLE PRECISION D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
263: COMPLEX*16 Q( LDQ, * ), WORK( * )
264: * ..
265: *
266: * =====================================================================
267: *
268: * .. Local Scalars ..
269: INTEGER COLTYP, CURR, I, IDLMDA, INDX,
270: $ INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
271: * ..
272: * .. External Subroutines ..
273: EXTERNAL DLAED9, DLAEDA, DLAMRG, XERBLA, ZLACRM, ZLAED8
274: * ..
275: * .. Intrinsic Functions ..
276: INTRINSIC MAX, MIN
277: * ..
278: * .. Executable Statements ..
279: *
280: * Test the input parameters.
281: *
282: INFO = 0
283: *
284: * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
285: * INFO = -1
286: * ELSE IF( N.LT.0 ) THEN
287: IF( N.LT.0 ) THEN
288: INFO = -1
289: ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
290: INFO = -2
291: ELSE IF( QSIZ.LT.N ) THEN
292: INFO = -3
293: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
294: INFO = -9
295: END IF
296: IF( INFO.NE.0 ) THEN
297: CALL XERBLA( 'ZLAED7', -INFO )
298: RETURN
299: END IF
300: *
301: * Quick return if possible
302: *
303: IF( N.EQ.0 )
304: $ RETURN
305: *
306: * The following values are for bookkeeping purposes only. They are
307: * integer pointers which indicate the portion of the workspace
308: * used by a particular array in DLAED2 and SLAED3.
309: *
310: IZ = 1
311: IDLMDA = IZ + N
312: IW = IDLMDA + N
313: IQ = IW + N
314: *
315: INDX = 1
316: INDXC = INDX + N
317: COLTYP = INDXC + N
318: INDXP = COLTYP + N
319: *
320: * Form the z-vector which consists of the last row of Q_1 and the
321: * first row of Q_2.
322: *
323: PTR = 1 + 2**TLVLS
324: DO 10 I = 1, CURLVL - 1
325: PTR = PTR + 2**( TLVLS-I )
326: 10 CONTINUE
327: CURR = PTR + CURPBM
328: CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
329: $ GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
330: $ RWORK( IZ+N ), INFO )
331: *
332: * When solving the final problem, we no longer need the stored data,
333: * so we will overwrite the data from this level onto the previously
334: * used storage space.
335: *
336: IF( CURLVL.EQ.TLVLS ) THEN
337: QPTR( CURR ) = 1
338: PRMPTR( CURR ) = 1
339: GIVPTR( CURR ) = 1
340: END IF
341: *
342: * Sort and Deflate eigenvalues.
343: *
344: CALL ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
345: $ RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
346: $ IWORK( INDXP ), IWORK( INDX ), INDXQ,
347: $ PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
348: $ GIVCOL( 1, GIVPTR( CURR ) ),
349: $ GIVNUM( 1, GIVPTR( CURR ) ), INFO )
350: PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
351: GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
352: *
353: * Solve Secular Equation.
354: *
355: IF( K.NE.0 ) THEN
356: CALL DLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
357: $ RWORK( IDLMDA ), RWORK( IW ),
358: $ QSTORE( QPTR( CURR ) ), K, INFO )
359: CALL ZLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
360: $ LDQ, RWORK( IQ ) )
361: QPTR( CURR+1 ) = QPTR( CURR ) + K**2
362: IF( INFO.NE.0 ) THEN
363: RETURN
364: END IF
365: *
366: * Prepare the INDXQ sorting premutation.
367: *
368: N1 = K
369: N2 = N - K
370: CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
371: ELSE
372: QPTR( CURR+1 ) = QPTR( CURR )
373: DO 20 I = 1, N
374: INDXQ( I ) = I
375: 20 CONTINUE
376: END IF
377: *
378: RETURN
379: *
380: * End of ZLAED7
381: *
382: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaed7.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack