Annotation of rpl/lapack/lapack/dlaed7.f, revision 1.3
1.1 bertrand 1: SUBROUTINE DLAED7( ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q,
2: $ LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR,
3: $ PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK,
4: $ INFO )
5: *
6: * -- LAPACK routine (version 3.2) --
7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
9: * November 2006
10: *
11: * .. Scalar Arguments ..
12: INTEGER CURLVL, CURPBM, CUTPNT, ICOMPQ, INFO, LDQ, N,
13: $ QSIZ, TLVLS
14: DOUBLE PRECISION RHO
15: * ..
16: * .. Array Arguments ..
17: INTEGER GIVCOL( 2, * ), GIVPTR( * ), INDXQ( * ),
18: $ IWORK( * ), PERM( * ), PRMPTR( * ), QPTR( * )
19: DOUBLE PRECISION D( * ), GIVNUM( 2, * ), Q( LDQ, * ),
20: $ QSTORE( * ), WORK( * )
21: * ..
22: *
23: * Purpose
24: * =======
25: *
26: * DLAED7 computes the updated eigensystem of a diagonal
27: * matrix after modification by a rank-one symmetric matrix. This
28: * routine is used only for the eigenproblem which requires all
29: * eigenvalues and optionally eigenvectors of a dense symmetric matrix
30: * that has been reduced to tridiagonal form. DLAED1 handles
31: * the case in which all eigenvalues and eigenvectors of a symmetric
32: * tridiagonal matrix are desired.
33: *
34: * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
35: *
36: * where Z = Q'u, u is a vector of length N with ones in the
37: * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
38: *
39: * The eigenvectors of the original matrix are stored in Q, and the
40: * eigenvalues are in D. The algorithm consists of three stages:
41: *
42: * The first stage consists of deflating the size of the problem
43: * when there are multiple eigenvalues or if there is a zero in
44: * the Z vector. For each such occurence the dimension of the
45: * secular equation problem is reduced by one. This stage is
46: * performed by the routine DLAED8.
47: *
48: * The second stage consists of calculating the updated
49: * eigenvalues. This is done by finding the roots of the secular
50: * equation via the routine DLAED4 (as called by DLAED9).
51: * This routine also calculates the eigenvectors of the current
52: * problem.
53: *
54: * The final stage consists of computing the updated eigenvectors
55: * directly using the updated eigenvalues. The eigenvectors for
56: * the current problem are multiplied with the eigenvectors from
57: * the overall problem.
58: *
59: * Arguments
60: * =========
61: *
62: * ICOMPQ (input) INTEGER
63: * = 0: Compute eigenvalues only.
64: * = 1: Compute eigenvectors of original dense symmetric matrix
65: * also. On entry, Q contains the orthogonal matrix used
66: * to reduce the original matrix to tridiagonal form.
67: *
68: * N (input) INTEGER
69: * The dimension of the symmetric tridiagonal matrix. N >= 0.
70: *
71: * QSIZ (input) INTEGER
72: * The dimension of the orthogonal matrix used to reduce
73: * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
74: *
75: * TLVLS (input) INTEGER
76: * The total number of merging levels in the overall divide and
77: * conquer tree.
78: *
79: * CURLVL (input) INTEGER
80: * The current level in the overall merge routine,
81: * 0 <= CURLVL <= TLVLS.
82: *
83: * CURPBM (input) INTEGER
84: * The current problem in the current level in the overall
85: * merge routine (counting from upper left to lower right).
86: *
87: * D (input/output) DOUBLE PRECISION array, dimension (N)
88: * On entry, the eigenvalues of the rank-1-perturbed matrix.
89: * On exit, the eigenvalues of the repaired matrix.
90: *
91: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
92: * On entry, the eigenvectors of the rank-1-perturbed matrix.
93: * On exit, the eigenvectors of the repaired tridiagonal matrix.
94: *
95: * LDQ (input) INTEGER
96: * The leading dimension of the array Q. LDQ >= max(1,N).
97: *
98: * INDXQ (output) INTEGER array, dimension (N)
99: * The permutation which will reintegrate the subproblem just
100: * solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
101: * will be in ascending order.
102: *
103: * RHO (input) DOUBLE PRECISION
104: * The subdiagonal element used to create the rank-1
105: * modification.
106: *
107: * CUTPNT (input) INTEGER
108: * Contains the location of the last eigenvalue in the leading
109: * sub-matrix. min(1,N) <= CUTPNT <= N.
110: *
111: * QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1)
112: * Stores eigenvectors of submatrices encountered during
113: * divide and conquer, packed together. QPTR points to
114: * beginning of the submatrices.
115: *
116: * QPTR (input/output) INTEGER array, dimension (N+2)
117: * List of indices pointing to beginning of submatrices stored
118: * in QSTORE. The submatrices are numbered starting at the
119: * bottom left of the divide and conquer tree, from left to
120: * right and bottom to top.
121: *
122: * PRMPTR (input) INTEGER array, dimension (N lg N)
123: * Contains a list of pointers which indicate where in PERM a
124: * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
125: * indicates the size of the permutation and also the size of
126: * the full, non-deflated problem.
127: *
128: * PERM (input) INTEGER array, dimension (N lg N)
129: * Contains the permutations (from deflation and sorting) to be
130: * applied to each eigenblock.
131: *
132: * GIVPTR (input) INTEGER array, dimension (N lg N)
133: * Contains a list of pointers which indicate where in GIVCOL a
134: * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
135: * indicates the number of Givens rotations.
136: *
137: * GIVCOL (input) INTEGER array, dimension (2, N lg N)
138: * Each pair of numbers indicates a pair of columns to take place
139: * in a Givens rotation.
140: *
141: * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
142: * Each number indicates the S value to be used in the
143: * corresponding Givens rotation.
144: *
145: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N)
146: *
147: * IWORK (workspace) INTEGER array, dimension (4*N)
148: *
149: * INFO (output) INTEGER
150: * = 0: successful exit.
151: * < 0: if INFO = -i, the i-th argument had an illegal value.
152: * > 0: if INFO = 1, an eigenvalue did not converge
153: *
154: * Further Details
155: * ===============
156: *
157: * Based on contributions by
158: * Jeff Rutter, Computer Science Division, University of California
159: * at Berkeley, USA
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION ONE, ZERO
165: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 )
166: * ..
167: * .. Local Scalars ..
168: INTEGER COLTYP, CURR, I, IDLMDA, INDX, INDXC, INDXP,
169: $ IQ2, IS, IW, IZ, K, LDQ2, N1, N2, PTR
170: * ..
171: * .. External Subroutines ..
172: EXTERNAL DGEMM, DLAED8, DLAED9, DLAEDA, DLAMRG, XERBLA
173: * ..
174: * .. Intrinsic Functions ..
175: INTRINSIC MAX, MIN
176: * ..
177: * .. Executable Statements ..
178: *
179: * Test the input parameters.
180: *
181: INFO = 0
182: *
183: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.1 ) THEN
184: INFO = -1
185: ELSE IF( N.LT.0 ) THEN
186: INFO = -2
187: ELSE IF( ICOMPQ.EQ.1 .AND. QSIZ.LT.N ) THEN
188: INFO = -4
189: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
190: INFO = -9
191: ELSE IF( MIN( 1, N ).GT.CUTPNT .OR. N.LT.CUTPNT ) THEN
192: INFO = -12
193: END IF
194: IF( INFO.NE.0 ) THEN
195: CALL XERBLA( 'DLAED7', -INFO )
196: RETURN
197: END IF
198: *
199: * Quick return if possible
200: *
201: IF( N.EQ.0 )
202: $ RETURN
203: *
204: * The following values are for bookkeeping purposes only. They are
205: * integer pointers which indicate the portion of the workspace
206: * used by a particular array in DLAED8 and DLAED9.
207: *
208: IF( ICOMPQ.EQ.1 ) THEN
209: LDQ2 = QSIZ
210: ELSE
211: LDQ2 = N
212: END IF
213: *
214: IZ = 1
215: IDLMDA = IZ + N
216: IW = IDLMDA + N
217: IQ2 = IW + N
218: IS = IQ2 + N*LDQ2
219: *
220: INDX = 1
221: INDXC = INDX + N
222: COLTYP = INDXC + N
223: INDXP = COLTYP + N
224: *
225: * Form the z-vector which consists of the last row of Q_1 and the
226: * first row of Q_2.
227: *
228: PTR = 1 + 2**TLVLS
229: DO 10 I = 1, CURLVL - 1
230: PTR = PTR + 2**( TLVLS-I )
231: 10 CONTINUE
232: CURR = PTR + CURPBM
233: CALL DLAEDA( N, TLVLS, CURLVL, CURPBM, PRMPTR, PERM, GIVPTR,
234: $ GIVCOL, GIVNUM, QSTORE, QPTR, WORK( IZ ),
235: $ WORK( IZ+N ), INFO )
236: *
237: * When solving the final problem, we no longer need the stored data,
238: * so we will overwrite the data from this level onto the previously
239: * used storage space.
240: *
241: IF( CURLVL.EQ.TLVLS ) THEN
242: QPTR( CURR ) = 1
243: PRMPTR( CURR ) = 1
244: GIVPTR( CURR ) = 1
245: END IF
246: *
247: * Sort and Deflate eigenvalues.
248: *
249: CALL DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, CUTPNT,
250: $ WORK( IZ ), WORK( IDLMDA ), WORK( IQ2 ), LDQ2,
251: $ WORK( IW ), PERM( PRMPTR( CURR ) ), GIVPTR( CURR+1 ),
252: $ GIVCOL( 1, GIVPTR( CURR ) ),
253: $ GIVNUM( 1, GIVPTR( CURR ) ), IWORK( INDXP ),
254: $ IWORK( INDX ), INFO )
255: PRMPTR( CURR+1 ) = PRMPTR( CURR ) + N
256: GIVPTR( CURR+1 ) = GIVPTR( CURR+1 ) + GIVPTR( CURR )
257: *
258: * Solve Secular Equation.
259: *
260: IF( K.NE.0 ) THEN
261: CALL DLAED9( K, 1, K, N, D, WORK( IS ), K, RHO, WORK( IDLMDA ),
262: $ WORK( IW ), QSTORE( QPTR( CURR ) ), K, INFO )
263: IF( INFO.NE.0 )
264: $ GO TO 30
265: IF( ICOMPQ.EQ.1 ) THEN
266: CALL DGEMM( 'N', 'N', QSIZ, K, K, ONE, WORK( IQ2 ), LDQ2,
267: $ QSTORE( QPTR( CURR ) ), K, ZERO, Q, LDQ )
268: END IF
269: QPTR( CURR+1 ) = QPTR( CURR ) + K**2
270: *
271: * Prepare the INDXQ sorting permutation.
272: *
273: N1 = K
274: N2 = N - K
275: CALL DLAMRG( N1, N2, D, 1, -1, INDXQ )
276: ELSE
277: QPTR( CURR+1 ) = QPTR( CURR )
278: DO 20 I = 1, N
279: INDXQ( I ) = I
280: 20 CONTINUE
281: END IF
282: *
283: 30 CONTINUE
284: RETURN
285: *
286: * End of DLAED7
287: *
288: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed7.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack