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