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