LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
zsteqr.f
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
9 *> Download ZSTEQR + dependencies
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
dlaev2
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
zsteqr
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zlasr
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
zlaset
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
dlascl
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
zswap
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
dlae2
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:102
dlasrt
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88