Annotation of rpl/lapack/lapack/dlaed0.f, revision 1.3
1.1 bertrand 1: SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
2: $ WORK, IWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
11: * ..
12: * .. Array Arguments ..
13: INTEGER IWORK( * )
14: DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
15: $ WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLAED0 computes all eigenvalues and corresponding eigenvectors of a
22: * symmetric tridiagonal matrix using the divide and conquer method.
23: *
24: * Arguments
25: * =========
26: *
27: * ICOMPQ (input) INTEGER
28: * = 0: Compute eigenvalues only.
29: * = 1: Compute eigenvectors of original dense symmetric matrix
30: * also. On entry, Q contains the orthogonal matrix used
31: * to reduce the original matrix to tridiagonal form.
32: * = 2: Compute eigenvalues and eigenvectors of tridiagonal
33: * matrix.
34: *
35: * QSIZ (input) INTEGER
36: * The dimension of the orthogonal matrix used to reduce
37: * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
38: *
39: * N (input) INTEGER
40: * The dimension of the symmetric tridiagonal matrix. N >= 0.
41: *
42: * D (input/output) DOUBLE PRECISION array, dimension (N)
43: * On entry, the main diagonal of the tridiagonal matrix.
44: * On exit, its eigenvalues.
45: *
46: * E (input) DOUBLE PRECISION array, dimension (N-1)
47: * The off-diagonal elements of the tridiagonal matrix.
48: * On exit, E has been destroyed.
49: *
50: * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
51: * On entry, Q must contain an N-by-N orthogonal matrix.
52: * If ICOMPQ = 0 Q is not referenced.
53: * If ICOMPQ = 1 On entry, Q is a subset of the columns of the
54: * orthogonal matrix used to reduce the full
55: * matrix to tridiagonal form corresponding to
56: * the subset of the full matrix which is being
57: * decomposed at this time.
58: * If ICOMPQ = 2 On entry, Q will be the identity matrix.
59: * On exit, Q contains the eigenvectors of the
60: * tridiagonal matrix.
61: *
62: * LDQ (input) INTEGER
63: * The leading dimension of the array Q. If eigenvectors are
64: * desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
65: *
66: * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
67: * Referenced only when ICOMPQ = 1. Used to store parts of
68: * the eigenvector matrix when the updating matrix multiplies
69: * take place.
70: *
71: * LDQS (input) INTEGER
72: * The leading dimension of the array QSTORE. If ICOMPQ = 1,
73: * then LDQS >= max(1,N). In any case, LDQS >= 1.
74: *
75: * WORK (workspace) DOUBLE PRECISION array,
76: * If ICOMPQ = 0 or 1, the dimension of WORK must be at least
77: * 1 + 3*N + 2*N*lg N + 2*N**2
78: * ( lg( N ) = smallest integer k
79: * such that 2^k >= N )
80: * If ICOMPQ = 2, the dimension of WORK must be at least
81: * 4*N + N**2.
82: *
83: * IWORK (workspace) INTEGER array,
84: * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
85: * 6 + 6*N + 5*N*lg N.
86: * ( lg( N ) = smallest integer k
87: * such that 2^k >= N )
88: * If ICOMPQ = 2, the dimension of IWORK must be at least
89: * 3 + 5*N.
90: *
91: * INFO (output) INTEGER
92: * = 0: successful exit.
93: * < 0: if INFO = -i, the i-th argument had an illegal value.
94: * > 0: The algorithm failed to compute an eigenvalue while
95: * working on the submatrix lying in rows and columns
96: * INFO/(N+1) through mod(INFO,N+1).
97: *
98: * Further Details
99: * ===============
100: *
101: * Based on contributions by
102: * Jeff Rutter, Computer Science Division, University of California
103: * at Berkeley, USA
104: *
105: * =====================================================================
106: *
107: * .. Parameters ..
108: DOUBLE PRECISION ZERO, ONE, TWO
109: PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
110: * ..
111: * .. Local Scalars ..
112: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
113: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
114: $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
115: $ SPM2, SUBMAT, SUBPBS, TLVLS
116: DOUBLE PRECISION TEMP
117: * ..
118: * .. External Subroutines ..
119: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
120: $ XERBLA
121: * ..
122: * .. External Functions ..
123: INTEGER ILAENV
124: EXTERNAL ILAENV
125: * ..
126: * .. Intrinsic Functions ..
127: INTRINSIC ABS, DBLE, INT, LOG, MAX
128: * ..
129: * .. Executable Statements ..
130: *
131: * Test the input parameters.
132: *
133: INFO = 0
134: *
135: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
136: INFO = -1
137: ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
138: INFO = -2
139: ELSE IF( N.LT.0 ) THEN
140: INFO = -3
141: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
142: INFO = -7
143: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
144: INFO = -9
145: END IF
146: IF( INFO.NE.0 ) THEN
147: CALL XERBLA( 'DLAED0', -INFO )
148: RETURN
149: END IF
150: *
151: * Quick return if possible
152: *
153: IF( N.EQ.0 )
154: $ RETURN
155: *
156: SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
157: *
158: * Determine the size and placement of the submatrices, and save in
159: * the leading elements of IWORK.
160: *
161: IWORK( 1 ) = N
162: SUBPBS = 1
163: TLVLS = 0
164: 10 CONTINUE
165: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
166: DO 20 J = SUBPBS, 1, -1
167: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
168: IWORK( 2*J-1 ) = IWORK( J ) / 2
169: 20 CONTINUE
170: TLVLS = TLVLS + 1
171: SUBPBS = 2*SUBPBS
172: GO TO 10
173: END IF
174: DO 30 J = 2, SUBPBS
175: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
176: 30 CONTINUE
177: *
178: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
179: * using rank-1 modifications (cuts).
180: *
181: SPM1 = SUBPBS - 1
182: DO 40 I = 1, SPM1
183: SUBMAT = IWORK( I ) + 1
184: SMM1 = SUBMAT - 1
185: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
186: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
187: 40 CONTINUE
188: *
189: INDXQ = 4*N + 3
190: IF( ICOMPQ.NE.2 ) THEN
191: *
192: * Set up workspaces for eigenvalues only/accumulate new vectors
193: * routine
194: *
195: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
196: LGN = INT( TEMP )
197: IF( 2**LGN.LT.N )
198: $ LGN = LGN + 1
199: IF( 2**LGN.LT.N )
200: $ LGN = LGN + 1
201: IPRMPT = INDXQ + N + 1
202: IPERM = IPRMPT + N*LGN
203: IQPTR = IPERM + N*LGN
204: IGIVPT = IQPTR + N + 2
205: IGIVCL = IGIVPT + N*LGN
206: *
207: IGIVNM = 1
208: IQ = IGIVNM + 2*N*LGN
209: IWREM = IQ + N**2 + 1
210: *
211: * Initialize pointers
212: *
213: DO 50 I = 0, SUBPBS
214: IWORK( IPRMPT+I ) = 1
215: IWORK( IGIVPT+I ) = 1
216: 50 CONTINUE
217: IWORK( IQPTR ) = 1
218: END IF
219: *
220: * Solve each submatrix eigenproblem at the bottom of the divide and
221: * conquer tree.
222: *
223: CURR = 0
224: DO 70 I = 0, SPM1
225: IF( I.EQ.0 ) THEN
226: SUBMAT = 1
227: MATSIZ = IWORK( 1 )
228: ELSE
229: SUBMAT = IWORK( I ) + 1
230: MATSIZ = IWORK( I+1 ) - IWORK( I )
231: END IF
232: IF( ICOMPQ.EQ.2 ) THEN
233: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
234: $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
235: IF( INFO.NE.0 )
236: $ GO TO 130
237: ELSE
238: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
239: $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
240: $ INFO )
241: IF( INFO.NE.0 )
242: $ GO TO 130
243: IF( ICOMPQ.EQ.1 ) THEN
244: CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
245: $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
246: $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
247: $ LDQS )
248: END IF
249: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
250: CURR = CURR + 1
251: END IF
252: K = 1
253: DO 60 J = SUBMAT, IWORK( I+1 )
254: IWORK( INDXQ+J ) = K
255: K = K + 1
256: 60 CONTINUE
257: 70 CONTINUE
258: *
259: * Successively merge eigensystems of adjacent submatrices
260: * into eigensystem for the corresponding larger matrix.
261: *
262: * while ( SUBPBS > 1 )
263: *
264: CURLVL = 1
265: 80 CONTINUE
266: IF( SUBPBS.GT.1 ) THEN
267: SPM2 = SUBPBS - 2
268: DO 90 I = 0, SPM2, 2
269: IF( I.EQ.0 ) THEN
270: SUBMAT = 1
271: MATSIZ = IWORK( 2 )
272: MSD2 = IWORK( 1 )
273: CURPRB = 0
274: ELSE
275: SUBMAT = IWORK( I ) + 1
276: MATSIZ = IWORK( I+2 ) - IWORK( I )
277: MSD2 = MATSIZ / 2
278: CURPRB = CURPRB + 1
279: END IF
280: *
281: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
282: * into an eigensystem of size MATSIZ.
283: * DLAED1 is used only for the full eigensystem of a tridiagonal
284: * matrix.
285: * DLAED7 handles the cases in which eigenvalues only or eigenvalues
286: * and eigenvectors of a full symmetric matrix (which was reduced to
287: * tridiagonal form) are desired.
288: *
289: IF( ICOMPQ.EQ.2 ) THEN
290: CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
291: $ LDQ, IWORK( INDXQ+SUBMAT ),
292: $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
293: $ IWORK( SUBPBS+1 ), INFO )
294: ELSE
295: CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
296: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
297: $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
298: $ MSD2, WORK( IQ ), IWORK( IQPTR ),
299: $ IWORK( IPRMPT ), IWORK( IPERM ),
300: $ IWORK( IGIVPT ), IWORK( IGIVCL ),
301: $ WORK( IGIVNM ), WORK( IWREM ),
302: $ IWORK( SUBPBS+1 ), INFO )
303: END IF
304: IF( INFO.NE.0 )
305: $ GO TO 130
306: IWORK( I / 2+1 ) = IWORK( I+2 )
307: 90 CONTINUE
308: SUBPBS = SUBPBS / 2
309: CURLVL = CURLVL + 1
310: GO TO 80
311: END IF
312: *
313: * end while
314: *
315: * Re-merge the eigenvalues/vectors which were deflated at the final
316: * merge step.
317: *
318: IF( ICOMPQ.EQ.1 ) THEN
319: DO 100 I = 1, N
320: J = IWORK( INDXQ+I )
321: WORK( I ) = D( J )
322: CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
323: 100 CONTINUE
324: CALL DCOPY( N, WORK, 1, D, 1 )
325: ELSE IF( ICOMPQ.EQ.2 ) THEN
326: DO 110 I = 1, N
327: J = IWORK( INDXQ+I )
328: WORK( I ) = D( J )
329: CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
330: 110 CONTINUE
331: CALL DCOPY( N, WORK, 1, D, 1 )
332: CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
333: ELSE
334: DO 120 I = 1, N
335: J = IWORK( INDXQ+I )
336: WORK( I ) = D( J )
337: 120 CONTINUE
338: CALL DCOPY( N, WORK, 1, D, 1 )
339: END IF
340: GO TO 140
341: *
342: 130 CONTINUE
343: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
344: *
345: 140 CONTINUE
346: RETURN
347: *
348: * End of DLAED0
349: *
350: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack