Annotation of rpl/lapack/lapack/dlaed0.f, revision 1.18
1.18 ! bertrand 1: *> \brief \b DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
1.8 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 9: *> Download DLAED0 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f">
1.8 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
22: * WORK, IWORK, INFO )
1.15 bertrand 23: *
1.8 bertrand 24: * .. Scalar Arguments ..
25: * INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
26: * ..
27: * .. Array Arguments ..
28: * INTEGER IWORK( * )
29: * DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
30: * $ WORK( * )
31: * ..
1.15 bertrand 32: *
1.8 bertrand 33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
40: *> symmetric tridiagonal matrix using the divide and conquer method.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] ICOMPQ
47: *> \verbatim
48: *> ICOMPQ is INTEGER
49: *> = 0: Compute eigenvalues only.
50: *> = 1: Compute eigenvectors of original dense symmetric matrix
51: *> also. On entry, Q contains the orthogonal matrix used
52: *> to reduce the original matrix to tridiagonal form.
53: *> = 2: Compute eigenvalues and eigenvectors of tridiagonal
54: *> matrix.
55: *> \endverbatim
56: *>
57: *> \param[in] QSIZ
58: *> \verbatim
59: *> QSIZ is INTEGER
60: *> The dimension of the orthogonal matrix used to reduce
61: *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
62: *> \endverbatim
63: *>
64: *> \param[in] N
65: *> \verbatim
66: *> N is INTEGER
67: *> The dimension of the symmetric tridiagonal matrix. N >= 0.
68: *> \endverbatim
69: *>
70: *> \param[in,out] D
71: *> \verbatim
72: *> D is DOUBLE PRECISION array, dimension (N)
73: *> On entry, the main diagonal of the tridiagonal matrix.
74: *> On exit, its eigenvalues.
75: *> \endverbatim
76: *>
77: *> \param[in] E
78: *> \verbatim
79: *> E is DOUBLE PRECISION array, dimension (N-1)
80: *> The off-diagonal elements of the tridiagonal matrix.
81: *> On exit, E has been destroyed.
82: *> \endverbatim
83: *>
84: *> \param[in,out] Q
85: *> \verbatim
86: *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
87: *> On entry, Q must contain an N-by-N orthogonal matrix.
88: *> If ICOMPQ = 0 Q is not referenced.
89: *> If ICOMPQ = 1 On entry, Q is a subset of the columns of the
90: *> orthogonal matrix used to reduce the full
91: *> matrix to tridiagonal form corresponding to
92: *> the subset of the full matrix which is being
93: *> decomposed at this time.
94: *> If ICOMPQ = 2 On entry, Q will be the identity matrix.
95: *> On exit, Q contains the eigenvectors of the
96: *> tridiagonal matrix.
97: *> \endverbatim
98: *>
99: *> \param[in] LDQ
100: *> \verbatim
101: *> LDQ is INTEGER
102: *> The leading dimension of the array Q. If eigenvectors are
103: *> desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
104: *> \endverbatim
105: *>
106: *> \param[out] QSTORE
107: *> \verbatim
108: *> QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
109: *> Referenced only when ICOMPQ = 1. Used to store parts of
110: *> the eigenvector matrix when the updating matrix multiplies
111: *> take place.
112: *> \endverbatim
113: *>
114: *> \param[in] LDQS
115: *> \verbatim
116: *> LDQS is INTEGER
117: *> The leading dimension of the array QSTORE. If ICOMPQ = 1,
118: *> then LDQS >= max(1,N). In any case, LDQS >= 1.
119: *> \endverbatim
120: *>
121: *> \param[out] WORK
122: *> \verbatim
123: *> WORK is DOUBLE PRECISION array,
124: *> If ICOMPQ = 0 or 1, the dimension of WORK must be at least
125: *> 1 + 3*N + 2*N*lg N + 3*N**2
126: *> ( lg( N ) = smallest integer k
127: *> such that 2^k >= N )
128: *> If ICOMPQ = 2, the dimension of WORK must be at least
129: *> 4*N + N**2.
130: *> \endverbatim
131: *>
132: *> \param[out] IWORK
133: *> \verbatim
134: *> IWORK is INTEGER array,
135: *> If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
136: *> 6 + 6*N + 5*N*lg N.
137: *> ( lg( N ) = smallest integer k
138: *> such that 2^k >= N )
139: *> If ICOMPQ = 2, the dimension of IWORK must be at least
140: *> 3 + 5*N.
141: *> \endverbatim
142: *>
143: *> \param[out] INFO
144: *> \verbatim
145: *> INFO is INTEGER
146: *> = 0: successful exit.
147: *> < 0: if INFO = -i, the i-th argument had an illegal value.
148: *> > 0: The algorithm failed to compute an eigenvalue while
149: *> working on the submatrix lying in rows and columns
150: *> INFO/(N+1) through mod(INFO,N+1).
151: *> \endverbatim
152: *
153: * Authors:
154: * ========
155: *
1.15 bertrand 156: *> \author Univ. of Tennessee
157: *> \author Univ. of California Berkeley
158: *> \author Univ. of Colorado Denver
159: *> \author NAG Ltd.
1.8 bertrand 160: *
161: *> \ingroup auxOTHERcomputational
162: *
163: *> \par Contributors:
164: * ==================
165: *>
166: *> Jeff Rutter, Computer Science Division, University of California
167: *> at Berkeley, USA
168: *
169: * =====================================================================
1.1 bertrand 170: SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
171: $ WORK, IWORK, INFO )
172: *
1.18 ! bertrand 173: * -- LAPACK computational routine --
1.1 bertrand 174: * -- LAPACK is a software package provided by Univ. of Tennessee, --
175: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176: *
177: * .. Scalar Arguments ..
178: INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
179: * ..
180: * .. Array Arguments ..
181: INTEGER IWORK( * )
182: DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
183: $ WORK( * )
184: * ..
185: *
186: * =====================================================================
187: *
188: * .. Parameters ..
189: DOUBLE PRECISION ZERO, ONE, TWO
190: PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
191: * ..
192: * .. Local Scalars ..
193: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
194: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
195: $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
196: $ SPM2, SUBMAT, SUBPBS, TLVLS
197: DOUBLE PRECISION TEMP
198: * ..
199: * .. External Subroutines ..
200: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
201: $ XERBLA
202: * ..
203: * .. External Functions ..
204: INTEGER ILAENV
205: EXTERNAL ILAENV
206: * ..
207: * .. Intrinsic Functions ..
208: INTRINSIC ABS, DBLE, INT, LOG, MAX
209: * ..
210: * .. Executable Statements ..
211: *
212: * Test the input parameters.
213: *
214: INFO = 0
215: *
216: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
217: INFO = -1
218: ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
219: INFO = -2
220: ELSE IF( N.LT.0 ) THEN
221: INFO = -3
222: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
223: INFO = -7
224: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
225: INFO = -9
226: END IF
227: IF( INFO.NE.0 ) THEN
228: CALL XERBLA( 'DLAED0', -INFO )
229: RETURN
230: END IF
231: *
232: * Quick return if possible
233: *
234: IF( N.EQ.0 )
235: $ RETURN
236: *
237: SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
238: *
239: * Determine the size and placement of the submatrices, and save in
240: * the leading elements of IWORK.
241: *
242: IWORK( 1 ) = N
243: SUBPBS = 1
244: TLVLS = 0
245: 10 CONTINUE
246: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
247: DO 20 J = SUBPBS, 1, -1
248: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
249: IWORK( 2*J-1 ) = IWORK( J ) / 2
250: 20 CONTINUE
251: TLVLS = TLVLS + 1
252: SUBPBS = 2*SUBPBS
253: GO TO 10
254: END IF
255: DO 30 J = 2, SUBPBS
256: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
257: 30 CONTINUE
258: *
259: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
260: * using rank-1 modifications (cuts).
261: *
262: SPM1 = SUBPBS - 1
263: DO 40 I = 1, SPM1
264: SUBMAT = IWORK( I ) + 1
265: SMM1 = SUBMAT - 1
266: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
267: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
268: 40 CONTINUE
269: *
270: INDXQ = 4*N + 3
271: IF( ICOMPQ.NE.2 ) THEN
272: *
273: * Set up workspaces for eigenvalues only/accumulate new vectors
274: * routine
275: *
276: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
277: LGN = INT( TEMP )
278: IF( 2**LGN.LT.N )
279: $ LGN = LGN + 1
280: IF( 2**LGN.LT.N )
281: $ LGN = LGN + 1
282: IPRMPT = INDXQ + N + 1
283: IPERM = IPRMPT + N*LGN
284: IQPTR = IPERM + N*LGN
285: IGIVPT = IQPTR + N + 2
286: IGIVCL = IGIVPT + N*LGN
287: *
288: IGIVNM = 1
289: IQ = IGIVNM + 2*N*LGN
290: IWREM = IQ + N**2 + 1
291: *
292: * Initialize pointers
293: *
294: DO 50 I = 0, SUBPBS
295: IWORK( IPRMPT+I ) = 1
296: IWORK( IGIVPT+I ) = 1
297: 50 CONTINUE
298: IWORK( IQPTR ) = 1
299: END IF
300: *
301: * Solve each submatrix eigenproblem at the bottom of the divide and
302: * conquer tree.
303: *
304: CURR = 0
305: DO 70 I = 0, SPM1
306: IF( I.EQ.0 ) THEN
307: SUBMAT = 1
308: MATSIZ = IWORK( 1 )
309: ELSE
310: SUBMAT = IWORK( I ) + 1
311: MATSIZ = IWORK( I+1 ) - IWORK( I )
312: END IF
313: IF( ICOMPQ.EQ.2 ) THEN
314: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
315: $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
316: IF( INFO.NE.0 )
317: $ GO TO 130
318: ELSE
319: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
320: $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
321: $ INFO )
322: IF( INFO.NE.0 )
323: $ GO TO 130
324: IF( ICOMPQ.EQ.1 ) THEN
325: CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
326: $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
327: $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
328: $ LDQS )
329: END IF
330: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
331: CURR = CURR + 1
332: END IF
333: K = 1
334: DO 60 J = SUBMAT, IWORK( I+1 )
335: IWORK( INDXQ+J ) = K
336: K = K + 1
337: 60 CONTINUE
338: 70 CONTINUE
339: *
340: * Successively merge eigensystems of adjacent submatrices
341: * into eigensystem for the corresponding larger matrix.
342: *
343: * while ( SUBPBS > 1 )
344: *
345: CURLVL = 1
346: 80 CONTINUE
347: IF( SUBPBS.GT.1 ) THEN
348: SPM2 = SUBPBS - 2
349: DO 90 I = 0, SPM2, 2
350: IF( I.EQ.0 ) THEN
351: SUBMAT = 1
352: MATSIZ = IWORK( 2 )
353: MSD2 = IWORK( 1 )
354: CURPRB = 0
355: ELSE
356: SUBMAT = IWORK( I ) + 1
357: MATSIZ = IWORK( I+2 ) - IWORK( I )
358: MSD2 = MATSIZ / 2
359: CURPRB = CURPRB + 1
360: END IF
361: *
362: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
363: * into an eigensystem of size MATSIZ.
364: * DLAED1 is used only for the full eigensystem of a tridiagonal
365: * matrix.
366: * DLAED7 handles the cases in which eigenvalues only or eigenvalues
367: * and eigenvectors of a full symmetric matrix (which was reduced to
368: * tridiagonal form) are desired.
369: *
370: IF( ICOMPQ.EQ.2 ) THEN
371: CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
372: $ LDQ, IWORK( INDXQ+SUBMAT ),
373: $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
374: $ IWORK( SUBPBS+1 ), INFO )
375: ELSE
376: CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
377: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
378: $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
379: $ MSD2, WORK( IQ ), IWORK( IQPTR ),
380: $ IWORK( IPRMPT ), IWORK( IPERM ),
381: $ IWORK( IGIVPT ), IWORK( IGIVCL ),
382: $ WORK( IGIVNM ), WORK( IWREM ),
383: $ IWORK( SUBPBS+1 ), INFO )
384: END IF
385: IF( INFO.NE.0 )
386: $ GO TO 130
387: IWORK( I / 2+1 ) = IWORK( I+2 )
388: 90 CONTINUE
389: SUBPBS = SUBPBS / 2
390: CURLVL = CURLVL + 1
391: GO TO 80
392: END IF
393: *
394: * end while
395: *
396: * Re-merge the eigenvalues/vectors which were deflated at the final
397: * merge step.
398: *
399: IF( ICOMPQ.EQ.1 ) THEN
400: DO 100 I = 1, N
401: J = IWORK( INDXQ+I )
402: WORK( I ) = D( J )
403: CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
404: 100 CONTINUE
405: CALL DCOPY( N, WORK, 1, D, 1 )
406: ELSE IF( ICOMPQ.EQ.2 ) THEN
407: DO 110 I = 1, N
408: J = IWORK( INDXQ+I )
409: WORK( I ) = D( J )
410: CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
411: 110 CONTINUE
412: CALL DCOPY( N, WORK, 1, D, 1 )
413: CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
414: ELSE
415: DO 120 I = 1, N
416: J = IWORK( INDXQ+I )
417: WORK( I ) = D( J )
418: 120 CONTINUE
419: CALL DCOPY( N, WORK, 1, D, 1 )
420: END IF
421: GO TO 140
422: *
423: 130 CONTINUE
424: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
425: *
426: 140 CONTINUE
427: RETURN
428: *
429: * End of DLAED0
430: *
431: 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