 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ zsteqr()

 subroutine zsteqr ( character COMPZ, integer N, double precision, dimension( * ) D, double precision, dimension( * ) E, complex*16, dimension( ldz, * ) Z, integer LDZ, double precision, dimension( * ) WORK, integer INFO )

ZSTEQR

Purpose:
``` ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a
symmetric tridiagonal matrix using the implicit QL or QR method.
The eigenvectors of a full or band complex Hermitian matrix can also
be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this
matrix to tridiagonal form.```
Parameters
 [in] COMPZ ``` COMPZ is CHARACTER*1 = 'N': Compute eigenvalues only. = 'V': Compute eigenvalues and eigenvectors of the original Hermitian matrix. On entry, Z must contain the unitary matrix used to reduce the original matrix to tridiagonal form. = 'I': Compute eigenvalues and eigenvectors of the tridiagonal matrix. Z is initialized to the identity matrix.``` [in] N ``` N is INTEGER The order of the matrix. N >= 0.``` [in,out] D ``` D is DOUBLE PRECISION array, dimension (N) On entry, the diagonal elements of the tridiagonal matrix. On exit, if INFO = 0, the eigenvalues in ascending order.``` [in,out] E ``` E is DOUBLE PRECISION array, dimension (N-1) On entry, the (n-1) subdiagonal elements of the tridiagonal matrix. On exit, E has been destroyed.``` [in,out] Z ``` Z is COMPLEX*16 array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the unitary matrix used in the reduction to tridiagonal form. On exit, if INFO = 0, then if COMPZ = 'V', Z contains the orthonormal eigenvectors of the original Hermitian matrix, and if COMPZ = 'I', Z contains the orthonormal eigenvectors of the symmetric tridiagonal matrix. If COMPZ = 'N', then Z is not referenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= 1, and if eigenvectors are desired, then LDZ >= max(1,N).``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (max(1,2*N-2)) If COMPZ = 'N', then WORK is not referenced.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: the algorithm has failed to find all the eigenvalues in a total of 30*N iterations; if INFO = i, then i elements of E have not converged to zero; on exit, D and E contain the elements of a symmetric tridiagonal matrix which is unitarily similar to the original matrix.```

Definition at line 131 of file zsteqr.f.

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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function dlanst(NORM, N, D, E)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: dlanst.f:100
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 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 dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:102
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dlasrt(ID, N, D, INFO)
DLASRT sorts numbers in increasing or decreasing order.
Definition: dlasrt.f:88
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
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
Here is the call graph for this function:
Here is the caller graph for this function: