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