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

Download ZSTEQR + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: