LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ csteqr()

subroutine csteqr ( character  COMPZ,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  INFO 
)

CSTEQR

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

Purpose:
 CSTEQR 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 CHETRD or CHPTRD or CHBTRD 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 REAL 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 REAL 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 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 REAL 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.
Date
December 2016

Definition at line 134 of file csteqr.f.

134 *
135 * -- LAPACK computational routine (version 3.7.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 * December 2016
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 *
subroutine slae2(A, B, C, RT1, RT2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: slae2.f:104
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:90
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: clasr.f:202
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine slaev2(A, B, C, RT1, RT2, CS1, SN1)
SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition: slaev2.f:122
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:83
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
real function slanst(NORM, N, D, E)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric tridiagonal matrix.
Definition: slanst.f:102
Here is the call graph for this function:
Here is the caller graph for this function: