LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zstedc.f
Go to the documentation of this file.
1*> \brief \b ZSTEDC
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZSTEDC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstedc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstedc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstedc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
22* LRWORK, IWORK, LIWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER COMPZ
26* INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
27* ..
28* .. Array Arguments ..
29* INTEGER IWORK( * )
30* DOUBLE PRECISION D( * ), E( * ), RWORK( * )
31* COMPLEX*16 WORK( * ), Z( LDZ, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> ZSTEDC computes all eigenvalues and, optionally, eigenvectors of a
41*> symmetric tridiagonal matrix using the divide and conquer method.
42*> The eigenvectors of a full or band complex Hermitian matrix can also
43*> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
44*> matrix to tridiagonal form.
45*>
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] COMPZ
52*> \verbatim
53*> COMPZ is CHARACTER*1
54*> = 'N': Compute eigenvalues only.
55*> = 'I': Compute eigenvectors of tridiagonal matrix also.
56*> = 'V': Compute eigenvectors of original Hermitian matrix
57*> also. On entry, Z contains the unitary matrix used
58*> to reduce the original matrix to tridiagonal form.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The dimension of the symmetric tridiagonal matrix. N >= 0.
65*> \endverbatim
66*>
67*> \param[in,out] D
68*> \verbatim
69*> D is DOUBLE PRECISION array, dimension (N)
70*> On entry, the diagonal elements of the tridiagonal matrix.
71*> On exit, if INFO = 0, the eigenvalues in ascending order.
72*> \endverbatim
73*>
74*> \param[in,out] E
75*> \verbatim
76*> E is DOUBLE PRECISION array, dimension (N-1)
77*> On entry, the subdiagonal elements of the tridiagonal matrix.
78*> On exit, E has been destroyed.
79*> \endverbatim
80*>
81*> \param[in,out] Z
82*> \verbatim
83*> Z is COMPLEX*16 array, dimension (LDZ,N)
84*> On entry, if COMPZ = 'V', then Z contains the unitary
85*> matrix used in the reduction to tridiagonal form.
86*> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
87*> orthonormal eigenvectors of the original Hermitian matrix,
88*> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
89*> of the symmetric tridiagonal matrix.
90*> If COMPZ = 'N', then Z is not referenced.
91*> \endverbatim
92*>
93*> \param[in] LDZ
94*> \verbatim
95*> LDZ is INTEGER
96*> The leading dimension of the array Z. LDZ >= 1.
97*> If eigenvectors are desired, then LDZ >= max(1,N).
98*> \endverbatim
99*>
100*> \param[out] WORK
101*> \verbatim
102*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
103*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
104*> \endverbatim
105*>
106*> \param[in] LWORK
107*> \verbatim
108*> LWORK is INTEGER
109*> The dimension of the array WORK.
110*> If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1.
111*> If COMPZ = 'V' and N > 1, LWORK must be at least N*N.
112*> Note that for COMPZ = 'V', then if N is less than or
113*> equal to the minimum divide size, usually 25, then LWORK need
114*> only be 1.
115*>
116*> If LWORK = -1, then a workspace query is assumed; the routine
117*> only calculates the optimal sizes of the WORK, RWORK and
118*> IWORK arrays, returns these values as the first entries of
119*> the WORK, RWORK and IWORK arrays, and no error message
120*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
121*> \endverbatim
122*>
123*> \param[out] RWORK
124*> \verbatim
125*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
126*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
127*> \endverbatim
128*>
129*> \param[in] LRWORK
130*> \verbatim
131*> LRWORK is INTEGER
132*> The dimension of the array RWORK.
133*> If COMPZ = 'N' or N <= 1, LRWORK must be at least 1.
134*> If COMPZ = 'V' and N > 1, LRWORK must be at least
135*> 1 + 3*N + 2*N*lg N + 4*N**2 ,
136*> where lg( N ) = smallest integer k such
137*> that 2**k >= N.
138*> If COMPZ = 'I' and N > 1, LRWORK must be at least
139*> 1 + 4*N + 2*N**2 .
140*> Note that for COMPZ = 'I' or 'V', then if N is less than or
141*> equal to the minimum divide size, usually 25, then LRWORK
142*> need only be max(1,2*(N-1)).
143*>
144*> If LRWORK = -1, then a workspace query is assumed; the
145*> routine only calculates the optimal sizes of the WORK, RWORK
146*> and IWORK arrays, returns these values as the first entries
147*> of the WORK, RWORK and IWORK arrays, and no error message
148*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
149*> \endverbatim
150*>
151*> \param[out] IWORK
152*> \verbatim
153*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
154*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
155*> \endverbatim
156*>
157*> \param[in] LIWORK
158*> \verbatim
159*> LIWORK is INTEGER
160*> The dimension of the array IWORK.
161*> If COMPZ = 'N' or N <= 1, LIWORK must be at least 1.
162*> If COMPZ = 'V' or N > 1, LIWORK must be at least
163*> 6 + 6*N + 5*N*lg N.
164*> If COMPZ = 'I' or N > 1, LIWORK must be at least
165*> 3 + 5*N .
166*> Note that for COMPZ = 'I' or 'V', then if N is less than or
167*> equal to the minimum divide size, usually 25, then LIWORK
168*> need only be 1.
169*>
170*> If LIWORK = -1, then a workspace query is assumed; the
171*> routine only calculates the optimal sizes of the WORK, RWORK
172*> and IWORK arrays, returns these values as the first entries
173*> of the WORK, RWORK and IWORK arrays, and no error message
174*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
175*> \endverbatim
176*>
177*> \param[out] INFO
178*> \verbatim
179*> INFO is INTEGER
180*> = 0: successful exit.
181*> < 0: if INFO = -i, the i-th argument had an illegal value.
182*> > 0: The algorithm failed to compute an eigenvalue while
183*> working on the submatrix lying in rows and columns
184*> INFO/(N+1) through mod(INFO,N+1).
185*> \endverbatim
186*
187* Authors:
188* ========
189*
190*> \author Univ. of Tennessee
191*> \author Univ. of California Berkeley
192*> \author Univ. of Colorado Denver
193*> \author NAG Ltd.
194*
195*> \ingroup stedc
196*
197*> \par Contributors:
198* ==================
199*>
200*> Jeff Rutter, Computer Science Division, University of California
201*> at Berkeley, USA
202*
203* =====================================================================
204 SUBROUTINE zstedc( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK,
205 $ LRWORK, IWORK, LIWORK, INFO )
206*
207* -- LAPACK computational routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 CHARACTER COMPZ
213 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
214* ..
215* .. Array Arguments ..
216 INTEGER IWORK( * )
217 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
218 COMPLEX*16 WORK( * ), Z( LDZ, * )
219* ..
220*
221* =====================================================================
222*
223* .. Parameters ..
224 DOUBLE PRECISION ZERO, ONE, TWO
225 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
226* ..
227* .. Local Scalars ..
228 LOGICAL LQUERY
229 INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN, LL,
230 $ lrwmin, lwmin, m, smlsiz, start
231 DOUBLE PRECISION EPS, ORGNRM, P, TINY
232* ..
233* .. External Functions ..
234 LOGICAL LSAME
235 INTEGER ILAENV
236 DOUBLE PRECISION DLAMCH, DLANST
237 EXTERNAL lsame, ilaenv, dlamch, dlanst
238* ..
239* .. External Subroutines ..
240 EXTERNAL dlascl, dlaset, dstedc, dsteqr, dsterf, xerbla,
242* ..
243* .. Intrinsic Functions ..
244 INTRINSIC abs, dble, int, log, max, mod, sqrt
245* ..
246* .. Executable Statements ..
247*
248* Test the input parameters.
249*
250 info = 0
251 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
252*
253 IF( lsame( compz, 'N' ) ) THEN
254 icompz = 0
255 ELSE IF( lsame( compz, 'V' ) ) THEN
256 icompz = 1
257 ELSE IF( lsame( compz, 'I' ) ) THEN
258 icompz = 2
259 ELSE
260 icompz = -1
261 END IF
262 IF( icompz.LT.0 ) THEN
263 info = -1
264 ELSE IF( n.LT.0 ) THEN
265 info = -2
266 ELSE IF( ( ldz.LT.1 ) .OR.
267 $ ( icompz.GT.0 .AND. ldz.LT.max( 1, n ) ) ) THEN
268 info = -6
269 END IF
270*
271 IF( info.EQ.0 ) THEN
272*
273* Compute the workspace requirements
274*
275 smlsiz = ilaenv( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 )
276 IF( n.LE.1 .OR. icompz.EQ.0 ) THEN
277 lwmin = 1
278 liwmin = 1
279 lrwmin = 1
280 ELSE IF( n.LE.smlsiz ) THEN
281 lwmin = 1
282 liwmin = 1
283 lrwmin = 2*( n - 1 )
284 ELSE IF( icompz.EQ.1 ) THEN
285 lgn = int( log( dble( n ) ) / log( two ) )
286 IF( 2**lgn.LT.n )
287 $ lgn = lgn + 1
288 IF( 2**lgn.LT.n )
289 $ lgn = lgn + 1
290 lwmin = n*n
291 lrwmin = 1 + 3*n + 2*n*lgn + 4*n**2
292 liwmin = 6 + 6*n + 5*n*lgn
293 ELSE IF( icompz.EQ.2 ) THEN
294 lwmin = 1
295 lrwmin = 1 + 4*n + 2*n**2
296 liwmin = 3 + 5*n
297 END IF
298 work( 1 ) = lwmin
299 rwork( 1 ) = lrwmin
300 iwork( 1 ) = liwmin
301*
302 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
303 info = -8
304 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
305 info = -10
306 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
307 info = -12
308 END IF
309 END IF
310*
311 IF( info.NE.0 ) THEN
312 CALL xerbla( 'ZSTEDC', -info )
313 RETURN
314 ELSE IF( lquery ) THEN
315 RETURN
316 END IF
317*
318* Quick return if possible
319*
320 IF( n.EQ.0 )
321 $ RETURN
322 IF( n.EQ.1 ) THEN
323 IF( icompz.NE.0 )
324 $ z( 1, 1 ) = one
325 RETURN
326 END IF
327*
328* If the following conditional clause is removed, then the routine
329* will use the Divide and Conquer routine to compute only the
330* eigenvalues, which requires (3N + 3N**2) real workspace and
331* (2 + 5N + 2N lg(N)) integer workspace.
332* Since on many architectures DSTERF is much faster than any other
333* algorithm for finding eigenvalues only, it is used here
334* as the default. If the conditional clause is removed, then
335* information on the size of workspace needs to be changed.
336*
337* If COMPZ = 'N', use DSTERF to compute the eigenvalues.
338*
339 IF( icompz.EQ.0 ) THEN
340 CALL dsterf( n, d, e, info )
341 GO TO 70
342 END IF
343*
344* If N is smaller than the minimum divide size (SMLSIZ+1), then
345* solve the problem with another solver.
346*
347 IF( n.LE.smlsiz ) THEN
348*
349 CALL zsteqr( compz, n, d, e, z, ldz, rwork, info )
350*
351 ELSE
352*
353* If COMPZ = 'I', we simply call DSTEDC instead.
354*
355 IF( icompz.EQ.2 ) THEN
356 CALL dlaset( 'Full', n, n, zero, one, rwork, n )
357 ll = n*n + 1
358 CALL dstedc( 'I', n, d, e, rwork, n,
359 $ rwork( ll ), lrwork-ll+1, iwork, liwork, info )
360 DO 20 j = 1, n
361 DO 10 i = 1, n
362 z( i, j ) = rwork( ( j-1 )*n+i )
363 10 CONTINUE
364 20 CONTINUE
365 GO TO 70
366 END IF
367*
368* From now on, only option left to be handled is COMPZ = 'V',
369* i.e. ICOMPZ = 1.
370*
371* Scale.
372*
373 orgnrm = dlanst( 'M', n, d, e )
374 IF( orgnrm.EQ.zero )
375 $ GO TO 70
376*
377 eps = dlamch( 'Epsilon' )
378*
379 start = 1
380*
381* while ( START <= N )
382*
383 30 CONTINUE
384 IF( start.LE.n ) THEN
385*
386* Let FINISH be the position of the next subdiagonal entry
387* such that E( FINISH ) <= TINY or FINISH = N if no such
388* subdiagonal exists. The matrix identified by the elements
389* between START and FINISH constitutes an independent
390* sub-problem.
391*
392 finish = start
393 40 CONTINUE
394 IF( finish.LT.n ) THEN
395 tiny = eps*sqrt( abs( d( finish ) ) )*
396 $ sqrt( abs( d( finish+1 ) ) )
397 IF( abs( e( finish ) ).GT.tiny ) THEN
398 finish = finish + 1
399 GO TO 40
400 END IF
401 END IF
402*
403* (Sub) Problem determined. Compute its size and solve it.
404*
405 m = finish - start + 1
406 IF( m.GT.smlsiz ) THEN
407*
408* Scale.
409*
410 orgnrm = dlanst( 'M', m, d( start ), e( start ) )
411 CALL dlascl( 'G', 0, 0, orgnrm, one, m, 1, d( start ), m,
412 $ info )
413 CALL dlascl( 'G', 0, 0, orgnrm, one, m-1, 1, e( start ),
414 $ m-1, info )
415*
416 CALL zlaed0( n, m, d( start ), e( start ), z( 1, start ),
417 $ ldz, work, n, rwork, iwork, info )
418 IF( info.GT.0 ) THEN
419 info = ( info / ( m+1 )+start-1 )*( n+1 ) +
420 $ mod( info, ( m+1 ) ) + start - 1
421 GO TO 70
422 END IF
423*
424* Scale back.
425*
426 CALL dlascl( 'G', 0, 0, one, orgnrm, m, 1, d( start ), m,
427 $ info )
428*
429 ELSE
430 CALL dsteqr( 'I', m, d( start ), e( start ), rwork, m,
431 $ rwork( m*m+1 ), info )
432 CALL zlacrm( n, m, z( 1, start ), ldz, rwork, m, work, n,
433 $ rwork( m*m+1 ) )
434 CALL zlacpy( 'A', n, m, work, n, z( 1, start ), ldz )
435 IF( info.GT.0 ) THEN
436 info = start*( n+1 ) + finish
437 GO TO 70
438 END IF
439 END IF
440*
441 start = finish + 1
442 GO TO 30
443 END IF
444*
445* endwhile
446*
447*
448* Use Selection Sort to minimize swaps of eigenvectors
449*
450 DO 60 ii = 2, n
451 i = ii - 1
452 k = i
453 p = d( i )
454 DO 50 j = ii, n
455 IF( d( j ).LT.p ) THEN
456 k = j
457 p = d( j )
458 END IF
459 50 CONTINUE
460 IF( k.NE.i ) THEN
461 d( k ) = d( i )
462 d( i ) = p
463 CALL zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
464 END IF
465 60 CONTINUE
466 END IF
467*
468 70 CONTINUE
469 work( 1 ) = lwmin
470 rwork( 1 ) = lrwmin
471 iwork( 1 ) = liwmin
472*
473 RETURN
474*
475* End of ZSTEDC
476*
477 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
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
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
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine zstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZSTEDC
Definition zstedc.f:206
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:182
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81