LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.

Definition at line 131 of file csteqr.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 REAL D( * ), E( * ), WORK( * )
143 COMPLEX Z( LDZ, * )
144* ..
145*
146* =====================================================================
147*
148* .. Parameters ..
149 REAL ZERO, ONE, TWO, THREE
150 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
151 $ three = 3.0e0 )
152 COMPLEX CZERO, CONE
153 parameter( czero = ( 0.0e0, 0.0e0 ),
154 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL 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 REAL SLAMCH, SLANST, SLAPY2
168 EXTERNAL lsame, slamch, slanst, slapy2
169* ..
170* .. External Subroutines ..
171 EXTERNAL claset, clasr, cswap, slae2, slaev2, slartg,
172 $ slascl, slasrt, xerbla
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( 'CSTEQR', -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 = slamch( 'E' )
219 eps2 = eps**2
220 safmin = slamch( '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 claset( '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 = slanst( '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 slascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
278 $ info )
279 CALL slascl( '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 slascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
284 $ info )
285 CALL slascl( '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 SLAE2 or SLAEV2
322* to compute its eigensystem.
323*
324 IF( m.EQ.l+1 ) THEN
325 IF( icompz.GT.0 ) THEN
326 CALL slaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
327 work( l ) = c
328 work( n-1+l ) = s
329 CALL clasr( 'R', 'V', 'B', n, 2, work( l ),
330 $ work( n-1+l ), z( 1, l ), ldz )
331 ELSE
332 CALL slae2( 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 = slapy2( 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 slartg( 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 clasr( '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 SLAE2 or SLAEV2
429* to compute its eigensystem.
430*
431 IF( m.EQ.l-1 ) THEN
432 IF( icompz.GT.0 ) THEN
433 CALL slaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
434 work( m ) = c
435 work( n-1+m ) = s
436 CALL clasr( 'R', 'V', 'F', n, 2, work( m ),
437 $ work( n-1+m ), z( 1, l-1 ), ldz )
438 ELSE
439 CALL slae2( 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 = slapy2( 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 slartg( 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 clasr( '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 slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
517 $ d( lsv ), n, info )
518 CALL slascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
519 $ n, info )
520 ELSE IF( iscale.EQ.2 ) THEN
521 CALL slascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
522 $ d( lsv ), n, info )
523 CALL slascl( '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 slasrt( '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 cswap( n, z( 1, i ), 1, z( 1, k ), 1 )
566 END IF
567 180 CONTINUE
568 END IF
569 RETURN
570*
571* End of CSTEQR
572*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slae2
subroutine slae2(a, b, c, rt1, rt2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition slae2.f:102
slaev2
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:120
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slanst
real function slanst(norm, n, d, e)
SLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slanst.f:100
slapy2
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
slascl
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:143
claset
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:106
clasr
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:200
slasrt
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:88
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: