 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dsteqr()

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

DSTEQR

Purpose:
``` DSTEQR 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 symmetric matrix can also be found
if DSYTRD or DSPTRD or DSBTRD 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 symmetric matrix. On entry, Z must contain the orthogonal 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 DOUBLE PRECISION array, dimension (LDZ, N) On entry, if COMPZ = 'V', then Z contains the orthogonal 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 symmetric 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 orthogonally similar to the original matrix.```

Definition at line 130 of file dsteqr.f.

131 *
132 * -- LAPACK computational routine --
133 * -- LAPACK is a software package provided by Univ. of Tennessee, --
134 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135 *
136 * .. Scalar Arguments ..
137  CHARACTER COMPZ
138  INTEGER INFO, LDZ, N
139 * ..
140 * .. Array Arguments ..
141  DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
142 * ..
143 *
144 * =====================================================================
145 *
146 * .. Parameters ..
147  DOUBLE PRECISION ZERO, ONE, TWO, THREE
148  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
149  \$ three = 3.0d0 )
150  INTEGER MAXIT
151  parameter( maxit = 30 )
152 * ..
153 * .. Local Scalars ..
154  INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
155  \$ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
156  \$ NM1, NMAXIT
157  DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
158  \$ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
159 * ..
160 * .. External Functions ..
161  LOGICAL LSAME
162  DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
163  EXTERNAL lsame, dlamch, dlanst, dlapy2
164 * ..
165 * .. External Subroutines ..
166  EXTERNAL dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr,
167  \$ dlasrt, dswap, xerbla
168 * ..
169 * .. Intrinsic Functions ..
170  INTRINSIC abs, max, sign, sqrt
171 * ..
172 * .. Executable Statements ..
173 *
174 * Test the input parameters.
175 *
176  info = 0
177 *
178  IF( lsame( compz, 'N' ) ) THEN
179  icompz = 0
180  ELSE IF( lsame( compz, 'V' ) ) THEN
181  icompz = 1
182  ELSE IF( lsame( compz, 'I' ) ) THEN
183  icompz = 2
184  ELSE
185  icompz = -1
186  END IF
187  IF( icompz.LT.0 ) THEN
188  info = -1
189  ELSE IF( n.LT.0 ) THEN
190  info = -2
191  ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
192  \$ n ) ) ) THEN
193  info = -6
194  END IF
195  IF( info.NE.0 ) THEN
196  CALL xerbla( 'DSTEQR', -info )
197  RETURN
198  END IF
199 *
200 * Quick return if possible
201 *
202  IF( n.EQ.0 )
203  \$ RETURN
204 *
205  IF( n.EQ.1 ) THEN
206  IF( icompz.EQ.2 )
207  \$ z( 1, 1 ) = one
208  RETURN
209  END IF
210 *
211 * Determine the unit roundoff and over/underflow thresholds.
212 *
213  eps = dlamch( 'E' )
214  eps2 = eps**2
215  safmin = dlamch( 'S' )
216  safmax = one / safmin
217  ssfmax = sqrt( safmax ) / three
218  ssfmin = sqrt( safmin ) / eps2
219 *
220 * Compute the eigenvalues and eigenvectors of the tridiagonal
221 * matrix.
222 *
223  IF( icompz.EQ.2 )
224  \$ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
225 *
226  nmaxit = n*maxit
227  jtot = 0
228 *
229 * Determine where the matrix splits and choose QL or QR iteration
230 * for each block, according to whether top or bottom diagonal
231 * element is smaller.
232 *
233  l1 = 1
234  nm1 = n - 1
235 *
236  10 CONTINUE
237  IF( l1.GT.n )
238  \$ GO TO 160
239  IF( l1.GT.1 )
240  \$ e( l1-1 ) = zero
241  IF( l1.LE.nm1 ) THEN
242  DO 20 m = l1, nm1
243  tst = abs( e( m ) )
244  IF( tst.EQ.zero )
245  \$ GO TO 30
246  IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
247  \$ 1 ) ) ) )*eps ) THEN
248  e( m ) = zero
249  GO TO 30
250  END IF
251  20 CONTINUE
252  END IF
253  m = n
254 *
255  30 CONTINUE
256  l = l1
257  lsv = l
258  lend = m
259  lendsv = lend
260  l1 = m + 1
261  IF( lend.EQ.l )
262  \$ GO TO 10
263 *
264 * Scale submatrix in rows and columns L to LEND
265 *
266  anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
267  iscale = 0
268  IF( anorm.EQ.zero )
269  \$ GO TO 10
270  IF( anorm.GT.ssfmax ) THEN
271  iscale = 1
272  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
273  \$ info )
274  CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
275  \$ info )
276  ELSE IF( anorm.LT.ssfmin ) THEN
277  iscale = 2
278  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
279  \$ info )
280  CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
281  \$ info )
282  END IF
283 *
284 * Choose between QL and QR iteration
285 *
286  IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
287  lend = lsv
288  l = lendsv
289  END IF
290 *
291  IF( lend.GT.l ) THEN
292 *
293 * QL Iteration
294 *
295 * Look for small subdiagonal element.
296 *
297  40 CONTINUE
298  IF( l.NE.lend ) THEN
299  lendm1 = lend - 1
300  DO 50 m = l, lendm1
301  tst = abs( e( m ) )**2
302  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
303  \$ safmin )GO TO 60
304  50 CONTINUE
305  END IF
306 *
307  m = lend
308 *
309  60 CONTINUE
310  IF( m.LT.lend )
311  \$ e( m ) = zero
312  p = d( l )
313  IF( m.EQ.l )
314  \$ GO TO 80
315 *
316 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
317 * to compute its eigensystem.
318 *
319  IF( m.EQ.l+1 ) THEN
320  IF( icompz.GT.0 ) THEN
321  CALL dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
322  work( l ) = c
323  work( n-1+l ) = s
324  CALL dlasr( 'R', 'V', 'B', n, 2, work( l ),
325  \$ work( n-1+l ), z( 1, l ), ldz )
326  ELSE
327  CALL dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
328  END IF
329  d( l ) = rt1
330  d( l+1 ) = rt2
331  e( l ) = zero
332  l = l + 2
333  IF( l.LE.lend )
334  \$ GO TO 40
335  GO TO 140
336  END IF
337 *
338  IF( jtot.EQ.nmaxit )
339  \$ GO TO 140
340  jtot = jtot + 1
341 *
342 * Form shift.
343 *
344  g = ( d( l+1 )-p ) / ( two*e( l ) )
345  r = dlapy2( g, one )
346  g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
347 *
348  s = one
349  c = one
350  p = zero
351 *
352 * Inner loop
353 *
354  mm1 = m - 1
355  DO 70 i = mm1, l, -1
356  f = s*e( i )
357  b = c*e( i )
358  CALL dlartg( g, f, c, s, r )
359  IF( i.NE.m-1 )
360  \$ e( i+1 ) = r
361  g = d( i+1 ) - p
362  r = ( d( i )-g )*s + two*c*b
363  p = s*r
364  d( i+1 ) = g + p
365  g = c*r - b
366 *
367 * If eigenvectors are desired, then save rotations.
368 *
369  IF( icompz.GT.0 ) THEN
370  work( i ) = c
371  work( n-1+i ) = -s
372  END IF
373 *
374  70 CONTINUE
375 *
376 * If eigenvectors are desired, then apply saved rotations.
377 *
378  IF( icompz.GT.0 ) THEN
379  mm = m - l + 1
380  CALL dlasr( 'R', 'V', 'B', n, mm, work( l ), work( n-1+l ),
381  \$ z( 1, l ), ldz )
382  END IF
383 *
384  d( l ) = d( l ) - p
385  e( l ) = g
386  GO TO 40
387 *
388 * Eigenvalue found.
389 *
390  80 CONTINUE
391  d( l ) = p
392 *
393  l = l + 1
394  IF( l.LE.lend )
395  \$ GO TO 40
396  GO TO 140
397 *
398  ELSE
399 *
400 * QR Iteration
401 *
402 * Look for small superdiagonal element.
403 *
404  90 CONTINUE
405  IF( l.NE.lend ) THEN
406  lendp1 = lend + 1
407  DO 100 m = l, lendp1, -1
408  tst = abs( e( m-1 ) )**2
409  IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
410  \$ safmin )GO TO 110
411  100 CONTINUE
412  END IF
413 *
414  m = lend
415 *
416  110 CONTINUE
417  IF( m.GT.lend )
418  \$ e( m-1 ) = zero
419  p = d( l )
420  IF( m.EQ.l )
421  \$ GO TO 130
422 *
423 * If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
424 * to compute its eigensystem.
425 *
426  IF( m.EQ.l-1 ) THEN
427  IF( icompz.GT.0 ) THEN
428  CALL dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
429  work( m ) = c
430  work( n-1+m ) = s
431  CALL dlasr( 'R', 'V', 'F', n, 2, work( m ),
432  \$ work( n-1+m ), z( 1, l-1 ), ldz )
433  ELSE
434  CALL dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
435  END IF
436  d( l-1 ) = rt1
437  d( l ) = rt2
438  e( l-1 ) = zero
439  l = l - 2
440  IF( l.GE.lend )
441  \$ GO TO 90
442  GO TO 140
443  END IF
444 *
445  IF( jtot.EQ.nmaxit )
446  \$ GO TO 140
447  jtot = jtot + 1
448 *
449 * Form shift.
450 *
451  g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
452  r = dlapy2( g, one )
453  g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
454 *
455  s = one
456  c = one
457  p = zero
458 *
459 * Inner loop
460 *
461  lm1 = l - 1
462  DO 120 i = m, lm1
463  f = s*e( i )
464  b = c*e( i )
465  CALL dlartg( g, f, c, s, r )
466  IF( i.NE.m )
467  \$ e( i-1 ) = r
468  g = d( i ) - p
469  r = ( d( i+1 )-g )*s + two*c*b
470  p = s*r
471  d( i ) = g + p
472  g = c*r - b
473 *
474 * If eigenvectors are desired, then save rotations.
475 *
476  IF( icompz.GT.0 ) THEN
477  work( i ) = c
478  work( n-1+i ) = s
479  END IF
480 *
481  120 CONTINUE
482 *
483 * If eigenvectors are desired, then apply saved rotations.
484 *
485  IF( icompz.GT.0 ) THEN
486  mm = l - m + 1
487  CALL dlasr( 'R', 'V', 'F', n, mm, work( m ), work( n-1+m ),
488  \$ z( 1, m ), ldz )
489  END IF
490 *
491  d( l ) = d( l ) - p
492  e( lm1 ) = g
493  GO TO 90
494 *
495 * Eigenvalue found.
496 *
497  130 CONTINUE
498  d( l ) = p
499 *
500  l = l - 1
501  IF( l.GE.lend )
502  \$ GO TO 90
503  GO TO 140
504 *
505  END IF
506 *
507 * Undo scaling if necessary
508 *
509  140 CONTINUE
510  IF( iscale.EQ.1 ) THEN
511  CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
512  \$ d( lsv ), n, info )
513  CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
514  \$ n, info )
515  ELSE IF( iscale.EQ.2 ) THEN
516  CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
517  \$ d( lsv ), n, info )
518  CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
519  \$ n, info )
520  END IF
521 *
522 * Check for no convergence to an eigenvalue after a total
523 * of N*MAXIT iterations.
524 *
525  IF( jtot.LT.nmaxit )
526  \$ GO TO 10
527  DO 150 i = 1, n - 1
528  IF( e( i ).NE.zero )
529  \$ info = info + 1
530  150 CONTINUE
531  GO TO 190
532 *
533 * Order eigenvalues and eigenvectors.
534 *
535  160 CONTINUE
536  IF( icompz.EQ.0 ) THEN
537 *
538 * Use Quick Sort
539 *
540  CALL dlasrt( 'I', n, d, info )
541 *
542  ELSE
543 *
544 * Use Selection Sort to minimize swaps of eigenvectors
545 *
546  DO 180 ii = 2, n
547  i = ii - 1
548  k = i
549  p = d( i )
550  DO 170 j = ii, n
551  IF( d( j ).LT.p ) THEN
552  k = j
553  p = d( j )
554  END IF
555  170 CONTINUE
556  IF( k.NE.i ) THEN
557  d( k ) = d( i )
558  d( i ) = p
559  CALL dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
560  END IF
561  180 CONTINUE
562  END IF
563 *
564  190 CONTINUE
565  RETURN
566 *
567 * End of DSTEQR
568 *
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 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 dlae2(A, B, C, RT1, RT2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: dlae2.f:102
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: dlasr.f:199
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 dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: