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