Annotation of rpl/lapack/lapack/dlaed0.f, revision 1.12
1.11 bertrand 1: *> \brief \b DLAED0 used by sstedc. 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: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
22: * WORK, IWORK, INFO )
23: *
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: * ..
32: *
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: *
156: *> \author Univ. of Tennessee
157: *> \author Univ. of California Berkeley
158: *> \author Univ. of Colorado Denver
159: *> \author NAG Ltd.
160: *
1.11 bertrand 161: *> \date September 2012
1.8 bertrand 162: *
163: *> \ingroup auxOTHERcomputational
164: *
165: *> \par Contributors:
166: * ==================
167: *>
168: *> Jeff Rutter, Computer Science Division, University of California
169: *> at Berkeley, USA
170: *
171: * =====================================================================
1.1 bertrand 172: SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
173: $ WORK, IWORK, INFO )
174: *
1.11 bertrand 175: * -- LAPACK computational routine (version 3.4.2) --
1.1 bertrand 176: * -- LAPACK is a software package provided by Univ. of Tennessee, --
177: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 178: * September 2012
1.1 bertrand 179: *
180: * .. Scalar Arguments ..
181: INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
182: * ..
183: * .. Array Arguments ..
184: INTEGER IWORK( * )
185: DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
186: $ WORK( * )
187: * ..
188: *
189: * =====================================================================
190: *
191: * .. Parameters ..
192: DOUBLE PRECISION ZERO, ONE, TWO
193: PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
194: * ..
195: * .. Local Scalars ..
196: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
197: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
198: $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
199: $ SPM2, SUBMAT, SUBPBS, TLVLS
200: DOUBLE PRECISION TEMP
201: * ..
202: * .. External Subroutines ..
203: EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
204: $ XERBLA
205: * ..
206: * .. External Functions ..
207: INTEGER ILAENV
208: EXTERNAL ILAENV
209: * ..
210: * .. Intrinsic Functions ..
211: INTRINSIC ABS, DBLE, INT, LOG, MAX
212: * ..
213: * .. Executable Statements ..
214: *
215: * Test the input parameters.
216: *
217: INFO = 0
218: *
219: IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
220: INFO = -1
221: ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
222: INFO = -2
223: ELSE IF( N.LT.0 ) THEN
224: INFO = -3
225: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
226: INFO = -7
227: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
228: INFO = -9
229: END IF
230: IF( INFO.NE.0 ) THEN
231: CALL XERBLA( 'DLAED0', -INFO )
232: RETURN
233: END IF
234: *
235: * Quick return if possible
236: *
237: IF( N.EQ.0 )
238: $ RETURN
239: *
240: SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
241: *
242: * Determine the size and placement of the submatrices, and save in
243: * the leading elements of IWORK.
244: *
245: IWORK( 1 ) = N
246: SUBPBS = 1
247: TLVLS = 0
248: 10 CONTINUE
249: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
250: DO 20 J = SUBPBS, 1, -1
251: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
252: IWORK( 2*J-1 ) = IWORK( J ) / 2
253: 20 CONTINUE
254: TLVLS = TLVLS + 1
255: SUBPBS = 2*SUBPBS
256: GO TO 10
257: END IF
258: DO 30 J = 2, SUBPBS
259: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
260: 30 CONTINUE
261: *
262: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
263: * using rank-1 modifications (cuts).
264: *
265: SPM1 = SUBPBS - 1
266: DO 40 I = 1, SPM1
267: SUBMAT = IWORK( I ) + 1
268: SMM1 = SUBMAT - 1
269: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
270: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
271: 40 CONTINUE
272: *
273: INDXQ = 4*N + 3
274: IF( ICOMPQ.NE.2 ) THEN
275: *
276: * Set up workspaces for eigenvalues only/accumulate new vectors
277: * routine
278: *
279: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
280: LGN = INT( TEMP )
281: IF( 2**LGN.LT.N )
282: $ LGN = LGN + 1
283: IF( 2**LGN.LT.N )
284: $ LGN = LGN + 1
285: IPRMPT = INDXQ + N + 1
286: IPERM = IPRMPT + N*LGN
287: IQPTR = IPERM + N*LGN
288: IGIVPT = IQPTR + N + 2
289: IGIVCL = IGIVPT + N*LGN
290: *
291: IGIVNM = 1
292: IQ = IGIVNM + 2*N*LGN
293: IWREM = IQ + N**2 + 1
294: *
295: * Initialize pointers
296: *
297: DO 50 I = 0, SUBPBS
298: IWORK( IPRMPT+I ) = 1
299: IWORK( IGIVPT+I ) = 1
300: 50 CONTINUE
301: IWORK( IQPTR ) = 1
302: END IF
303: *
304: * Solve each submatrix eigenproblem at the bottom of the divide and
305: * conquer tree.
306: *
307: CURR = 0
308: DO 70 I = 0, SPM1
309: IF( I.EQ.0 ) THEN
310: SUBMAT = 1
311: MATSIZ = IWORK( 1 )
312: ELSE
313: SUBMAT = IWORK( I ) + 1
314: MATSIZ = IWORK( I+1 ) - IWORK( I )
315: END IF
316: IF( ICOMPQ.EQ.2 ) THEN
317: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
318: $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
319: IF( INFO.NE.0 )
320: $ GO TO 130
321: ELSE
322: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
323: $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
324: $ INFO )
325: IF( INFO.NE.0 )
326: $ GO TO 130
327: IF( ICOMPQ.EQ.1 ) THEN
328: CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
329: $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
330: $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
331: $ LDQS )
332: END IF
333: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
334: CURR = CURR + 1
335: END IF
336: K = 1
337: DO 60 J = SUBMAT, IWORK( I+1 )
338: IWORK( INDXQ+J ) = K
339: K = K + 1
340: 60 CONTINUE
341: 70 CONTINUE
342: *
343: * Successively merge eigensystems of adjacent submatrices
344: * into eigensystem for the corresponding larger matrix.
345: *
346: * while ( SUBPBS > 1 )
347: *
348: CURLVL = 1
349: 80 CONTINUE
350: IF( SUBPBS.GT.1 ) THEN
351: SPM2 = SUBPBS - 2
352: DO 90 I = 0, SPM2, 2
353: IF( I.EQ.0 ) THEN
354: SUBMAT = 1
355: MATSIZ = IWORK( 2 )
356: MSD2 = IWORK( 1 )
357: CURPRB = 0
358: ELSE
359: SUBMAT = IWORK( I ) + 1
360: MATSIZ = IWORK( I+2 ) - IWORK( I )
361: MSD2 = MATSIZ / 2
362: CURPRB = CURPRB + 1
363: END IF
364: *
365: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
366: * into an eigensystem of size MATSIZ.
367: * DLAED1 is used only for the full eigensystem of a tridiagonal
368: * matrix.
369: * DLAED7 handles the cases in which eigenvalues only or eigenvalues
370: * and eigenvectors of a full symmetric matrix (which was reduced to
371: * tridiagonal form) are desired.
372: *
373: IF( ICOMPQ.EQ.2 ) THEN
374: CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
375: $ LDQ, IWORK( INDXQ+SUBMAT ),
376: $ E( SUBMAT+MSD2-1 ), MSD2, WORK,
377: $ IWORK( SUBPBS+1 ), INFO )
378: ELSE
379: CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
380: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
381: $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
382: $ MSD2, WORK( IQ ), IWORK( IQPTR ),
383: $ IWORK( IPRMPT ), IWORK( IPERM ),
384: $ IWORK( IGIVPT ), IWORK( IGIVCL ),
385: $ WORK( IGIVNM ), WORK( IWREM ),
386: $ IWORK( SUBPBS+1 ), INFO )
387: END IF
388: IF( INFO.NE.0 )
389: $ GO TO 130
390: IWORK( I / 2+1 ) = IWORK( I+2 )
391: 90 CONTINUE
392: SUBPBS = SUBPBS / 2
393: CURLVL = CURLVL + 1
394: GO TO 80
395: END IF
396: *
397: * end while
398: *
399: * Re-merge the eigenvalues/vectors which were deflated at the final
400: * merge step.
401: *
402: IF( ICOMPQ.EQ.1 ) THEN
403: DO 100 I = 1, N
404: J = IWORK( INDXQ+I )
405: WORK( I ) = D( J )
406: CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
407: 100 CONTINUE
408: CALL DCOPY( N, WORK, 1, D, 1 )
409: ELSE IF( ICOMPQ.EQ.2 ) THEN
410: DO 110 I = 1, N
411: J = IWORK( INDXQ+I )
412: WORK( I ) = D( J )
413: CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
414: 110 CONTINUE
415: CALL DCOPY( N, WORK, 1, D, 1 )
416: CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
417: ELSE
418: DO 120 I = 1, N
419: J = IWORK( INDXQ+I )
420: WORK( I ) = D( J )
421: 120 CONTINUE
422: CALL DCOPY( N, WORK, 1, D, 1 )
423: END IF
424: GO TO 140
425: *
426: 130 CONTINUE
427: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
428: *
429: 140 CONTINUE
430: RETURN
431: *
432: * End of DLAED0
433: *
434: 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