Annotation of rpl/lapack/lapack/zlaed7.f, revision 1.16
1.12 bertrand 1: *> \brief \b ZLAED7 used by sstedc. 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: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
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 )
25: *
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: * ..
37: *
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: *
237: *> \author Univ. of Tennessee
238: *> \author Univ. of California Berkeley
239: *> \author Univ. of Colorado Denver
240: *> \author NAG Ltd.
241: *
1.15 bertrand 242: *> \date June 2016
1.9 bertrand 243: *
244: *> \ingroup complex16OTHERcomputational
245: *
246: * =====================================================================
1.1 bertrand 247: SUBROUTINE ZLAED7( N, CUTPNT, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
248: $ LDQ, RHO, INDXQ, QSTORE, QPTR, PRMPTR, PERM,
249: $ GIVPTR, GIVCOL, GIVNUM, WORK, RWORK, IWORK,
250: $ INFO )
251: *
1.15 bertrand 252: * -- LAPACK computational routine (version 3.6.1) --
1.1 bertrand 253: * -- LAPACK is a software package provided by Univ. of Tennessee, --
254: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 bertrand 255: * June 2016
1.1 bertrand 256: *
257: * .. Scalar Arguments ..
258: INTEGER CURLVL, CURPBM, CUTPNT, INFO, LDQ, N, QSIZ,
259: $ TLVLS
260: DOUBLE PRECISION RHO
261: * ..
262: * .. Array Arguments ..
263: INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
264: $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
265: DOUBLE PRECISION D( * ), GIVNUM( 2, * ), QSTORE( * ), RWORK( * )
266: COMPLEX*16 Q( LDQ, * ), WORK( * )
267: * ..
268: *
269: * =====================================================================
270: *
271: * .. Local Scalars ..
272: INTEGER COLTYP, CURR, I, IDLMDA, INDX,
273: $ INDXC, INDXP, IQ, IW, IZ, K, N1, N2, PTR
274: * ..
275: * .. External Subroutines ..
276: EXTERNAL DLAED9, DLAEDA, DLAMRG, XERBLA, ZLACRM, ZLAED8
277: * ..
278: * .. Intrinsic Functions ..
279: INTRINSIC MAX, MIN
280: * ..
281: * .. Executable Statements ..
282: *
283: * Test the input parameters.
284: *
285: INFO = 0
286: *
287: * IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
288: * INFO = -1
289: * ELSE IF( N.LT.0 ) THEN
290: IF( N.LT.0 ) THEN
291: INFO = -1
292: ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
293: INFO = -2
294: ELSE IF( QSIZ.LT.N ) THEN
295: INFO = -3
296: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
297: INFO = -9
298: END IF
299: IF( INFO.NE.0 ) THEN
300: CALL XERBLA( 'ZLAED7', -INFO )
301: RETURN
302: END IF
303: *
304: * Quick return if possible
305: *
306: IF( N.EQ.0 )
307: $ RETURN
308: *
309: * The following values are for bookkeeping purposes only. They are
310: * integer pointers which indicate the portion of the workspace
311: * used by a particular array in DLAED2 and SLAED3.
312: *
313: IZ = 1
314: IDLMDA = IZ + N
315: IW = IDLMDA + N
316: IQ = IW + N
317: *
318: INDX = 1
319: INDXC = INDX + N
320: COLTYP = INDXC + N
321: INDXP = COLTYP + N
322: *
323: * Form the z-vector which consists of the last row of Q_1 and the
324: * first row of Q_2.
325: *
326: PTR = 1 + 2**TLVLS
327: DO 10 I = 1, CURLVL - 1
328: PTR = PTR + 2**( TLVLS-I )
329: 10 CONTINUE
330: CURR = PTR + CURPBM
331: CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
332: $ GIVCOL, GIVNUM, QSTORE, QPTR, RWORK( IZ ),
333: $ RWORK( IZ+N ), INFO )
334: *
335: * When solving the final problem, we no longer need the stored data,
336: * so we will overwrite the data from this level onto the previously
337: * used storage space.
338: *
339: IF( CURLVL.EQ.TLVLS ) THEN
340: QPTR( CURR ) = 1
341: PRMPTR( CURR ) = 1
342: GIVPTR( CURR ) = 1
343: END IF
344: *
345: * Sort and Deflate eigenvalues.
346: *
347: CALL ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, RWORK( IZ ),
348: $ RWORK( IDLMDA ), WORK, QSIZ, RWORK( IW ),
349: $ IWORK( INDXP ), IWORK( INDX ), INDXQ,
350: $ PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
351: $ GIVCOL( 1, GIVPTR( CURR ) ),
352: $ GIVNUM( 1, GIVPTR( CURR ) ), INFO )
353: PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
354: GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
355: *
356: * Solve Secular Equation.
357: *
358: IF( K.NE.0 ) THEN
359: CALL DLAED9( K, 1, K, N, D, RWORK( IQ ), K, RHO,
360: $ RWORK( IDLMDA ), RWORK( IW ),
361: $ QSTORE( QPTR( CURR ) ), K, INFO )
362: CALL ZLACRM( QSIZ, K, WORK, QSIZ, QSTORE( QPTR( CURR ) ), K, Q,
363: $ LDQ, RWORK( IQ ) )
364: QPTR( CURR+1 ) = QPTR( CURR ) + K**2
365: IF( INFO.NE.0 ) THEN
366: RETURN
367: END IF
368: *
369: * Prepare the INDXQ sorting premutation.
370: *
371: N1 = K
372: N2 = N - K
373: CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
374: ELSE
375: QPTR( CURR+1 ) = QPTR( CURR )
376: DO 20 I = 1, N
377: INDXQ( I ) = I
378: 20 CONTINUE
379: END IF
380: *
381: RETURN
382: *
383: * End of ZLAED7
384: *
385: 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