Annotation of rpl/lapack/lapack/zlaed0.f, revision 1.18
1.18 ! bertrand 1: *> \brief \b ZLAED0 used by ZSTEDC. 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 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">
1.8 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
22: * IWORK, INFO )
1.15 bertrand 23: *
1.8 bertrand 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: * ..
1.15 bertrand 32: *
1.8 bertrand 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: *
1.15 bertrand 135: *> \author Univ. of Tennessee
136: *> \author Univ. of California Berkeley
137: *> \author Univ. of Colorado Denver
138: *> \author NAG Ltd.
1.8 bertrand 139: *
140: *> \ingroup complex16OTHERcomputational
141: *
142: * =====================================================================
1.1 bertrand 143: SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
144: $ IWORK, INFO )
145: *
1.18 ! bertrand 146: * -- LAPACK computational routine --
1.1 bertrand 147: * -- LAPACK is a software package provided by Univ. of Tennessee, --
148: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149: *
150: * .. Scalar Arguments ..
151: INTEGER INFO, LDQ, LDQS, N, QSIZ
152: * ..
153: * .. Array Arguments ..
154: INTEGER IWORK( * )
155: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
156: COMPLEX*16 Q( LDQ, * ), QSTORE( LDQS, * )
157: * ..
158: *
159: * =====================================================================
160: *
161: * Warning: N could be as big as QSIZ!
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION TWO
165: PARAMETER ( TWO = 2.D+0 )
166: * ..
167: * .. Local Scalars ..
168: INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
169: $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
170: $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
171: $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
172: DOUBLE PRECISION TEMP
173: * ..
174: * .. External Subroutines ..
175: EXTERNAL DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
176: * ..
177: * .. External Functions ..
178: INTEGER ILAENV
179: EXTERNAL ILAENV
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC ABS, DBLE, INT, LOG, MAX
183: * ..
184: * .. Executable Statements ..
185: *
186: * Test the input parameters.
187: *
188: INFO = 0
189: *
190: * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
191: * INFO = -1
192: * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
193: * $ THEN
194: IF( QSIZ.LT.MAX( 0, N ) ) THEN
195: INFO = -1
196: ELSE IF( N.LT.0 ) THEN
197: INFO = -2
198: ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
199: INFO = -6
200: ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
201: INFO = -8
202: END IF
203: IF( INFO.NE.0 ) THEN
204: CALL XERBLA( 'ZLAED0', -INFO )
205: RETURN
206: END IF
207: *
208: * Quick return if possible
209: *
210: IF( N.EQ.0 )
211: $ RETURN
212: *
213: SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
214: *
215: * Determine the size and placement of the submatrices, and save in
216: * the leading elements of IWORK.
217: *
218: IWORK( 1 ) = N
219: SUBPBS = 1
220: TLVLS = 0
221: 10 CONTINUE
222: IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
223: DO 20 J = SUBPBS, 1, -1
224: IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
225: IWORK( 2*J-1 ) = IWORK( J ) / 2
226: 20 CONTINUE
227: TLVLS = TLVLS + 1
228: SUBPBS = 2*SUBPBS
229: GO TO 10
230: END IF
231: DO 30 J = 2, SUBPBS
232: IWORK( J ) = IWORK( J ) + IWORK( J-1 )
233: 30 CONTINUE
234: *
235: * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
236: * using rank-1 modifications (cuts).
237: *
238: SPM1 = SUBPBS - 1
239: DO 40 I = 1, SPM1
240: SUBMAT = IWORK( I ) + 1
241: SMM1 = SUBMAT - 1
242: D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
243: D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
244: 40 CONTINUE
245: *
246: INDXQ = 4*N + 3
247: *
248: * Set up workspaces for eigenvalues only/accumulate new vectors
249: * routine
250: *
251: TEMP = LOG( DBLE( N ) ) / LOG( TWO )
252: LGN = INT( TEMP )
253: IF( 2**LGN.LT.N )
254: $ LGN = LGN + 1
255: IF( 2**LGN.LT.N )
256: $ LGN = LGN + 1
257: IPRMPT = INDXQ + N + 1
258: IPERM = IPRMPT + N*LGN
259: IQPTR = IPERM + N*LGN
260: IGIVPT = IQPTR + N + 2
261: IGIVCL = IGIVPT + N*LGN
262: *
263: IGIVNM = 1
264: IQ = IGIVNM + 2*N*LGN
265: IWREM = IQ + N**2 + 1
266: * Initialize pointers
267: DO 50 I = 0, SUBPBS
268: IWORK( IPRMPT+I ) = 1
269: IWORK( IGIVPT+I ) = 1
270: 50 CONTINUE
271: IWORK( IQPTR ) = 1
272: *
273: * Solve each submatrix eigenproblem at the bottom of the divide and
274: * conquer tree.
275: *
276: CURR = 0
277: DO 70 I = 0, SPM1
278: IF( I.EQ.0 ) THEN
279: SUBMAT = 1
280: MATSIZ = IWORK( 1 )
281: ELSE
282: SUBMAT = IWORK( I ) + 1
283: MATSIZ = IWORK( I+1 ) - IWORK( I )
284: END IF
285: LL = IQ - 1 + IWORK( IQPTR+CURR )
286: CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
287: $ RWORK( LL ), MATSIZ, RWORK, INFO )
288: CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
289: $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
290: $ RWORK( IWREM ) )
291: IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
292: CURR = CURR + 1
293: IF( INFO.GT.0 ) THEN
294: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
295: RETURN
296: END IF
297: K = 1
298: DO 60 J = SUBMAT, IWORK( I+1 )
299: IWORK( INDXQ+J ) = K
300: K = K + 1
301: 60 CONTINUE
302: 70 CONTINUE
303: *
304: * Successively merge eigensystems of adjacent submatrices
305: * into eigensystem for the corresponding larger matrix.
306: *
307: * while ( SUBPBS > 1 )
308: *
309: CURLVL = 1
310: 80 CONTINUE
311: IF( SUBPBS.GT.1 ) THEN
312: SPM2 = SUBPBS - 2
313: DO 90 I = 0, SPM2, 2
314: IF( I.EQ.0 ) THEN
315: SUBMAT = 1
316: MATSIZ = IWORK( 2 )
317: MSD2 = IWORK( 1 )
318: CURPRB = 0
319: ELSE
320: SUBMAT = IWORK( I ) + 1
321: MATSIZ = IWORK( I+2 ) - IWORK( I )
322: MSD2 = MATSIZ / 2
323: CURPRB = CURPRB + 1
324: END IF
325: *
326: * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
327: * into an eigensystem of size MATSIZ. ZLAED7 handles the case
328: * when the eigenvectors of a full or band Hermitian matrix (which
329: * was reduced to tridiagonal form) are desired.
330: *
331: * I am free to use Q as a valuable working space until Loop 150.
332: *
333: CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
334: $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
335: $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
336: $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
337: $ IWORK( IPERM ), IWORK( IGIVPT ),
338: $ IWORK( IGIVCL ), RWORK( IGIVNM ),
339: $ Q( 1, SUBMAT ), RWORK( IWREM ),
340: $ IWORK( SUBPBS+1 ), INFO )
341: IF( INFO.GT.0 ) THEN
342: INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
343: RETURN
344: END IF
345: IWORK( I / 2+1 ) = IWORK( I+2 )
346: 90 CONTINUE
347: SUBPBS = SUBPBS / 2
348: CURLVL = CURLVL + 1
349: GO TO 80
350: END IF
351: *
352: * end while
353: *
354: * Re-merge the eigenvalues/vectors which were deflated at the final
355: * merge step.
356: *
357: DO 100 I = 1, N
358: J = IWORK( INDXQ+I )
359: RWORK( I ) = D( J )
360: CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
361: 100 CONTINUE
362: CALL DCOPY( N, RWORK, 1, D, 1 )
363: *
364: RETURN
365: *
366: * End of ZLAED0
367: *
368: 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