LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zsteqr.f
Go to the documentation of this file.
1*> \brief \b ZSTEQR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zsteqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zsteqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zsteqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER COMPZ
25* INTEGER INFO, LDZ, N
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION D( * ), E( * ), WORK( * )
29* COMPLEX*16 Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a
39*> symmetric tridiagonal matrix using the implicit QL or QR method.
40*> The eigenvectors of a full or band complex Hermitian matrix can also
41*> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
42*> matrix to tridiagonal form.
43*> \endverbatim
44*
45* Arguments:
46* ==========
47*
48*> \param[in] COMPZ
49*> \verbatim
50*> COMPZ is CHARACTER*1
51*> = 'N': Compute eigenvalues only.
52*> = 'V': Compute eigenvalues and eigenvectors of the original
53*> Hermitian matrix. On entry, Z must contain the
54*> unitary matrix used to reduce the original matrix
55*> to tridiagonal form.
56*> = 'I': Compute eigenvalues and eigenvectors of the
57*> tridiagonal matrix. Z is initialized to the identity
58*> matrix.
59*> \endverbatim
60*>
61*> \param[in] N
62*> \verbatim
63*> N is INTEGER
64*> The order of the 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 (n-1) subdiagonal elements of the tridiagonal
78*> matrix.
79*> On exit, E has been destroyed.
80*> \endverbatim
81*>
82*> \param[in,out] Z
83*> \verbatim
84*> Z is COMPLEX*16 array, dimension (LDZ, N)
85*> On entry, if COMPZ = 'V', then Z contains the unitary
86*> matrix used in the reduction to tridiagonal form.
87*> On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
88*> orthonormal eigenvectors of the original Hermitian matrix,
89*> and if COMPZ = 'I', Z contains the orthonormal eigenvectors
90*> of the symmetric tridiagonal matrix.
91*> If COMPZ = 'N', then Z is not referenced.
92*> \endverbatim
93*>
94*> \param[in] LDZ
95*> \verbatim
96*> LDZ is INTEGER
97*> The leading dimension of the array Z. LDZ >= 1, and if
98*> eigenvectors are desired, then LDZ >= max(1,N).
99*> \endverbatim
100*>
101*> \param[out] WORK
102*> \verbatim
103*> WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2))
104*> If COMPZ = 'N', then WORK is not referenced.
105*> \endverbatim
106*>
107*> \param[out] INFO
108*> \verbatim
109*> INFO is INTEGER
110*> = 0: successful exit
111*> < 0: if INFO = -i, the i-th argument had an illegal value
112*> > 0: the algorithm has failed to find all the eigenvalues in
113*> a total of 30*N iterations; if INFO = i, then i
114*> elements of E have not converged to zero; on exit, D
115*> and E contain the elements of a symmetric tridiagonal
116*> matrix which is unitarily similar to the original
117*> matrix.
118*> \endverbatim
119*
120* Authors:
121* ========
122*
123*> \author Univ. of Tennessee
124*> \author Univ. of California Berkeley
125*> \author Univ. of Colorado Denver
126*> \author NAG Ltd.
127*
128*> \ingroup steqr
129*
130* =====================================================================
131 SUBROUTINE zsteqr( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
132*
133* -- LAPACK computational routine --
134* -- LAPACK is a software package provided by Univ. of Tennessee, --
135* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137* .. Scalar Arguments ..
138 CHARACTER COMPZ
139 INTEGER INFO, LDZ, N
140* ..
141* .. Array Arguments ..
142 DOUBLE PRECISION D( * ), E( * ), WORK( * )
143 COMPLEX*16 Z( LDZ, * )
144* ..
145*
146* =====================================================================
147*
148* .. Parameters ..
149 DOUBLE PRECISION ZERO, ONE, TWO, THREE
150 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
151 \$ three = 3.0d0 )
152 COMPLEX*16 CZERO, CONE
153 parameter( czero = ( 0.0d0, 0.0d0 ),
154 \$ cone = ( 1.0d0, 0.0d0 ) )
155 INTEGER MAXIT
156 parameter( maxit = 30 )
157* ..
158* .. Local Scalars ..
159 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
160 \$ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
161 \$ NM1, NMAXIT
162 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
163 \$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
164* ..
165* .. External Functions ..
166 LOGICAL LSAME
167 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
168 EXTERNAL lsame, dlamch, dlanst, dlapy2
169* ..
170* .. External Subroutines ..
171 EXTERNAL dlae2, dlaev2, dlartg, dlascl, dlasrt, xerbla,
172 \$ zlaset, zlasr, zswap
173* ..
174* .. Intrinsic Functions ..
175 INTRINSIC abs, max, sign, sqrt
176* ..
177* .. Executable Statements ..
178*
179* Test the input parameters.
180*
181 info = 0
182*
183 IF( lsame( compz, 'N' ) ) THEN
184 icompz = 0
185 ELSE IF( lsame( compz, 'V' ) ) THEN
186 icompz = 1
187 ELSE IF( lsame( compz, 'I' ) ) THEN
188 icompz = 2
189 ELSE
190 icompz = -1
191 END IF
192 IF( icompz.LT.0 ) THEN
193 info = -1
194 ELSE IF( n.LT.0 ) THEN
195 info = -2
196 ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
197 \$ n ) ) ) THEN
198 info = -6
199 END IF
200 IF( info.NE.0 ) THEN
201 CALL xerbla( 'ZSTEQR', -info )
202 RETURN
203 END IF
204*
205* Quick return if possible
206*
207 IF( n.EQ.0 )
208 \$ RETURN
209*
210 IF( n.EQ.1 ) THEN
211 IF( icompz.EQ.2 )
212 \$ z( 1, 1 ) = cone
213 RETURN
214 END IF
215*
216* Determine the unit roundoff and over/underflow thresholds.
217*
218 eps = dlamch( 'E' )
219 eps2 = eps**2
220 safmin = dlamch( 'S' )
221 safmax = one / safmin
222 ssfmax = sqrt( safmax ) / three
223 ssfmin = sqrt( safmin ) / eps2
224*
225* Compute the eigenvalues and eigenvectors of the tridiagonal
226* matrix.
227*
228 IF( icompz.EQ.2 )
229 \$ CALL zlaset( 'Full', n, n, czero, cone, z, ldz )
230*
231 nmaxit = n*maxit
232 jtot = 0
233*
234* Determine where the matrix splits and choose QL or QR iteration
235* for each block, according to whether top or bottom diagonal
236* element is smaller.
237*
238 l1 = 1
239 nm1 = n - 1
240*
241 10 CONTINUE
242 IF( l1.GT.n )
243 \$ GO TO 160
244 IF( l1.GT.1 )
245 \$ e( l1-1 ) = zero
246 IF( l1.LE.nm1 ) THEN
247 DO 20 m = l1, nm1
248 tst = abs( e( m ) )
249 IF( tst.EQ.zero )
250 \$ GO TO 30
251 IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
252 \$ 1 ) ) ) )*eps ) THEN
253 e( m ) = zero
254 GO TO 30
255 END IF
256 20 CONTINUE
257 END IF
258 m = n
259*
260 30 CONTINUE
261 l = l1
262 lsv = l
263 lend = m
264 lendsv = lend
265 l1 = m + 1
266 IF( lend.EQ.l )
267 \$ GO TO 10
268*
269* Scale submatrix in rows and columns L to LEND
270*
271 anorm = dlanst( 'I', lend-l+1, d( l ), e( l ) )
272 iscale = 0
273 IF( anorm.EQ.zero )
274 \$ GO TO 10
275 IF( anorm.GT.ssfmax ) THEN
276 iscale = 1
277 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
278 \$ info )
279 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
280 \$ info )
281 ELSE IF( anorm.LT.ssfmin ) THEN
282 iscale = 2
283 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
284 \$ info )
285 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
286 \$ info )
287 END IF
288*
289* Choose between QL and QR iteration
290*
291 IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
292 lend = lsv
293 l = lendsv
294 END IF
295*
296 IF( lend.GT.l ) THEN
297*
298* QL Iteration
299*
300* Look for small subdiagonal element.
301*
302 40 CONTINUE
303 IF( l.NE.lend ) THEN
304 lendm1 = lend - 1
305 DO 50 m = l, lendm1
306 tst = abs( e( m ) )**2
307 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
308 \$ safmin )GO TO 60
309 50 CONTINUE
310 END IF
311*
312 m = lend
313*
314 60 CONTINUE
315 IF( m.LT.lend )
316 \$ e( m ) = zero
317 p = d( l )
318 IF( m.EQ.l )
319 \$ GO TO 80
320*
321* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
322* to compute its eigensystem.
323*
324 IF( m.EQ.l+1 ) THEN
325 IF( icompz.GT.0 ) THEN
326 CALL dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
327 work( l ) = c
328 work( n-1+l ) = s
329 CALL zlasr( 'R', 'V', 'B', n, 2, work( l ),
330 \$ work( n-1+l ), z( 1, l ), ldz )
331 ELSE
332 CALL dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
333 END IF
334 d( l ) = rt1
335 d( l+1 ) = rt2
336 e( l ) = zero
337 l = l + 2
338 IF( l.LE.lend )
339 \$ GO TO 40
340 GO TO 140
341 END IF
342*
343 IF( jtot.EQ.nmaxit )
344 \$ GO TO 140
345 jtot = jtot + 1
346*
347* Form shift.
348*
349 g = ( d( l+1 )-p ) / ( two*e( l ) )
350 r = dlapy2( g, one )
351 g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
352*
353 s = one
354 c = one
355 p = zero
356*
357* Inner loop
358*
359 mm1 = m - 1
360 DO 70 i = mm1, l, -1
361 f = s*e( i )
362 b = c*e( i )
363 CALL dlartg( g, f, c, s, r )
364 IF( i.NE.m-1 )
365 \$ e( i+1 ) = r
366 g = d( i+1 ) - p
367 r = ( d( i )-g )*s + two*c*b
368 p = s*r
369 d( i+1 ) = g + p
370 g = c*r - b
371*
372* If eigenvectors are desired, then save rotations.
373*
374 IF( icompz.GT.0 ) THEN
375 work( i ) = c
376 work( n-1+i ) = -s
377 END IF
378*
379 70 CONTINUE
380*
381* If eigenvectors are desired, then apply saved rotations.
382*
383 IF( icompz.GT.0 ) THEN
384 mm = m - l + 1
385 CALL zlasr( 'R', 'V', 'B', n, mm, work( l ), work( n-1+l ),
386 \$ z( 1, l ), ldz )
387 END IF
388*
389 d( l ) = d( l ) - p
390 e( l ) = g
391 GO TO 40
392*
393* Eigenvalue found.
394*
395 80 CONTINUE
396 d( l ) = p
397*
398 l = l + 1
399 IF( l.LE.lend )
400 \$ GO TO 40
401 GO TO 140
402*
403 ELSE
404*
405* QR Iteration
406*
407* Look for small superdiagonal element.
408*
409 90 CONTINUE
410 IF( l.NE.lend ) THEN
411 lendp1 = lend + 1
412 DO 100 m = l, lendp1, -1
413 tst = abs( e( m-1 ) )**2
414 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
415 \$ safmin )GO TO 110
416 100 CONTINUE
417 END IF
418*
419 m = lend
420*
421 110 CONTINUE
422 IF( m.GT.lend )
423 \$ e( m-1 ) = zero
424 p = d( l )
425 IF( m.EQ.l )
426 \$ GO TO 130
427*
428* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
429* to compute its eigensystem.
430*
431 IF( m.EQ.l-1 ) THEN
432 IF( icompz.GT.0 ) THEN
433 CALL dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
434 work( m ) = c
435 work( n-1+m ) = s
436 CALL zlasr( 'R', 'V', 'F', n, 2, work( m ),
437 \$ work( n-1+m ), z( 1, l-1 ), ldz )
438 ELSE
439 CALL dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
440 END IF
441 d( l-1 ) = rt1
442 d( l ) = rt2
443 e( l-1 ) = zero
444 l = l - 2
445 IF( l.GE.lend )
446 \$ GO TO 90
447 GO TO 140
448 END IF
449*
450 IF( jtot.EQ.nmaxit )
451 \$ GO TO 140
452 jtot = jtot + 1
453*
454* Form shift.
455*
456 g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
457 r = dlapy2( g, one )
458 g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
459*
460 s = one
461 c = one
462 p = zero
463*
464* Inner loop
465*
466 lm1 = l - 1
467 DO 120 i = m, lm1
468 f = s*e( i )
469 b = c*e( i )
470 CALL dlartg( g, f, c, s, r )
471 IF( i.NE.m )
472 \$ e( i-1 ) = r
473 g = d( i ) - p
474 r = ( d( i+1 )-g )*s + two*c*b
475 p = s*r
476 d( i ) = g + p
477 g = c*r - b
478*
479* If eigenvectors are desired, then save rotations.
480*
481 IF( icompz.GT.0 ) THEN
482 work( i ) = c
483 work( n-1+i ) = s
484 END IF
485*
486 120 CONTINUE
487*
488* If eigenvectors are desired, then apply saved rotations.
489*
490 IF( icompz.GT.0 ) THEN
491 mm = l - m + 1
492 CALL zlasr( 'R', 'V', 'F', n, mm, work( m ), work( n-1+m ),
493 \$ z( 1, m ), ldz )
494 END IF
495*
496 d( l ) = d( l ) - p
497 e( lm1 ) = g
498 GO TO 90
499*
500* Eigenvalue found.
501*
502 130 CONTINUE
503 d( l ) = p
504*
505 l = l - 1
506 IF( l.GE.lend )
507 \$ GO TO 90
508 GO TO 140
509*
510 END IF
511*
512* Undo scaling if necessary
513*
514 140 CONTINUE
515 IF( iscale.EQ.1 ) THEN
516 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
517 \$ d( lsv ), n, info )
518 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
519 \$ n, info )
520 ELSE IF( iscale.EQ.2 ) THEN
521 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
522 \$ d( lsv ), n, info )
523 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
524 \$ n, info )
525 END IF
526*
527* Check for no convergence to an eigenvalue after a total
528* of N*MAXIT iterations.
529*
530 IF( jtot.EQ.nmaxit ) THEN
531 DO 150 i = 1, n - 1
532 IF( e( i ).NE.zero )
533 \$ info = info + 1
534 150 CONTINUE
535 RETURN
536 END IF
537 GO TO 10
538*
539* Order eigenvalues and eigenvectors.
540*
541 160 CONTINUE
542 IF( icompz.EQ.0 ) THEN
543*
544* Use Quick Sort
545*
546 CALL dlasrt( 'I', n, d, info )
547*
548 ELSE
549*
550* Use Selection Sort to minimize swaps of eigenvectors
551*
552 DO 180 ii = 2, n
553 i = ii - 1
554 k = i
555 p = d( i )
556 DO 170 j = ii, n
557 IF( d( j ).LT.p ) THEN
558 k = j
559 p = d( j )
560 END IF
561 170 CONTINUE
562 IF( k.NE.i ) THEN
563 d( k ) = d( i )
564 d( i ) = p
565 CALL zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
566 END IF
567 180 CONTINUE
568 END IF
569 RETURN
570*
571* End of ZSTEQR
572*
573 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
subroutine dlaev2(a, b, c, rt1, rt2, cs1, sn1)
DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition dlaev2.f:120
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zlasr(side, pivot, direct, m, n, c, s, a, lda)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition zlasr.f:200
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81