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

◆ zsteqr()

subroutine zsteqr ( character  compz,
integer  n,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
complex*16, dimension( ldz, * )  z,
integer  ldz,
double precision, dimension( * )  work,
integer  info 
)

ZSTEQR

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

Purpose:
 ZSTEQR 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 ZHETRD or ZHPTRD or ZHBTRD 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 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 COMPLEX*16 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 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 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 zsteqr.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 DOUBLE PRECISION D( * ), E( * ), WORK( * )
143 COMPLEX*16 Z( LDZ, * )
144* ..
145*
146* =====================================================================
147*
148* .. Parameters ..
149 DOUBLE PRECISION ZERO, ONE, TWO, THREE
150 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
151 $ three = 3.0d0 )
152 COMPLEX*16 CZERO, CONE
153 parameter( czero = ( 0.0d0, 0.0d0 ),
154 $ cone = ( 1.0d0, 0.0d0 ) )
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 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
168 EXTERNAL lsame, dlamch, dlanst, dlapy2
169* ..
170* .. External Subroutines ..
171 EXTERNAL dlae2, dlaev2, dlartg, dlascl, dlasrt, xerbla,
172 $ zlaset, zlasr, zswap
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( 'ZSTEQR', -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 = dlamch( 'E' )
219 eps2 = eps**2
220 safmin = dlamch( '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 zlaset( '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 = dlanst( '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 dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
278 $ info )
279 CALL dlascl( '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 dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
284 $ info )
285 CALL dlascl( '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 DLAE2 or SLAEV2
322* to compute its eigensystem.
323*
324 IF( m.EQ.l+1 ) THEN
325 IF( icompz.GT.0 ) THEN
326 CALL dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
327 work( l ) = c
328 work( n-1+l ) = s
329 CALL zlasr( 'R', 'V', 'B', n, 2, work( l ),
330 $ work( n-1+l ), z( 1, l ), ldz )
331 ELSE
332 CALL dlae2( 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 = dlapy2( 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 dlartg( 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 zlasr( '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 DLAE2 or SLAEV2
429* to compute its eigensystem.
430*
431 IF( m.EQ.l-1 ) THEN
432 IF( icompz.GT.0 ) THEN
433 CALL dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
434 work( m ) = c
435 work( n-1+m ) = s
436 CALL zlasr( 'R', 'V', 'F', n, 2, work( m ),
437 $ work( n-1+m ), z( 1, l-1 ), ldz )
438 ELSE
439 CALL dlae2( 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 = dlapy2( 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 dlartg( 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 zlasr( '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 dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
517 $ d( lsv ), n, info )
518 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
519 $ n, info )
520 ELSE IF( iscale.EQ.2 ) THEN
521 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
522 $ d( lsv ), n, info )
523 CALL dlascl( '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 dlasrt( '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 zswap( n, z( 1, i ), 1, z( 1, k ), 1 )
566 END IF
567 180 CONTINUE
568 END IF
569 RETURN
570*
571* End of ZSTEQR
572*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
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
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
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zlasr(side, pivot, direct, m, n, c, s, a, lda)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition zlasr.f:200
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: