LAPACK 3.12.0
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*> \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*> \ingroup laed0
141*
142* =====================================================================
143 SUBROUTINE zlaed0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
144 $ IWORK, INFO )
145*
146* -- LAPACK computational routine --
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
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:114
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:145
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:249
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131