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

◆ dsteqr()

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

DSTEQR

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

Purpose:
 DSTEQR 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 symmetric matrix can also be found
 if DSYTRD or DSPTRD or DSBTRD 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
                  symmetric matrix.  On entry, Z must contain the
                  orthogonal 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 DOUBLE PRECISION array, dimension (LDZ, N)
          On entry, if  COMPZ = 'V', then Z contains the orthogonal
          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 symmetric 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 orthogonally similar to the original
                matrix.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 130 of file dsteqr.f.

131*
132* -- LAPACK computational routine --
133* -- LAPACK is a software package provided by Univ. of Tennessee, --
134* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
135*
136* .. Scalar Arguments ..
137 CHARACTER COMPZ
138 INTEGER INFO, LDZ, N
139* ..
140* .. Array Arguments ..
141 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
142* ..
143*
144* =====================================================================
145*
146* .. Parameters ..
147 DOUBLE PRECISION ZERO, ONE, TWO, THREE
148 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
149 $ three = 3.0d0 )
150 INTEGER MAXIT
151 parameter( maxit = 30 )
152* ..
153* .. Local Scalars ..
154 INTEGER I, ICOMPZ, II, ISCALE, J, JTOT, K, L, L1, LEND,
155 $ LENDM1, LENDP1, LENDSV, LM1, LSV, M, MM, MM1,
156 $ NM1, NMAXIT
157 DOUBLE PRECISION ANORM, B, C, EPS, EPS2, F, G, P, R, RT1, RT2,
158 $ S, SAFMAX, SAFMIN, SSFMAX, SSFMIN, TST
159* ..
160* .. External Functions ..
161 LOGICAL LSAME
162 DOUBLE PRECISION DLAMCH, DLANST, DLAPY2
163 EXTERNAL lsame, dlamch, dlanst, dlapy2
164* ..
165* .. External Subroutines ..
166 EXTERNAL dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr,
167 $ dlasrt, dswap, xerbla
168* ..
169* .. Intrinsic Functions ..
170 INTRINSIC abs, max, sign, sqrt
171* ..
172* .. Executable Statements ..
173*
174* Test the input parameters.
175*
176 info = 0
177*
178 IF( lsame( compz, 'N' ) ) THEN
179 icompz = 0
180 ELSE IF( lsame( compz, 'V' ) ) THEN
181 icompz = 1
182 ELSE IF( lsame( compz, 'I' ) ) THEN
183 icompz = 2
184 ELSE
185 icompz = -1
186 END IF
187 IF( icompz.LT.0 ) THEN
188 info = -1
189 ELSE IF( n.LT.0 ) THEN
190 info = -2
191 ELSE IF( ( ldz.LT.1 ) .OR. ( icompz.GT.0 .AND. ldz.LT.max( 1,
192 $ n ) ) ) THEN
193 info = -6
194 END IF
195 IF( info.NE.0 ) THEN
196 CALL xerbla( 'DSTEQR', -info )
197 RETURN
198 END IF
199*
200* Quick return if possible
201*
202 IF( n.EQ.0 )
203 $ RETURN
204*
205 IF( n.EQ.1 ) THEN
206 IF( icompz.EQ.2 )
207 $ z( 1, 1 ) = one
208 RETURN
209 END IF
210*
211* Determine the unit roundoff and over/underflow thresholds.
212*
213 eps = dlamch( 'E' )
214 eps2 = eps**2
215 safmin = dlamch( 'S' )
216 safmax = one / safmin
217 ssfmax = sqrt( safmax ) / three
218 ssfmin = sqrt( safmin ) / eps2
219*
220* Compute the eigenvalues and eigenvectors of the tridiagonal
221* matrix.
222*
223 IF( icompz.EQ.2 )
224 $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
225*
226 nmaxit = n*maxit
227 jtot = 0
228*
229* Determine where the matrix splits and choose QL or QR iteration
230* for each block, according to whether top or bottom diagonal
231* element is smaller.
232*
233 l1 = 1
234 nm1 = n - 1
235*
236 10 CONTINUE
237 IF( l1.GT.n )
238 $ GO TO 160
239 IF( l1.GT.1 )
240 $ e( l1-1 ) = zero
241 IF( l1.LE.nm1 ) THEN
242 DO 20 m = l1, nm1
243 tst = abs( e( m ) )
244 IF( tst.EQ.zero )
245 $ GO TO 30
246 IF( tst.LE.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+
247 $ 1 ) ) ) )*eps ) THEN
248 e( m ) = zero
249 GO TO 30
250 END IF
251 20 CONTINUE
252 END IF
253 m = n
254*
255 30 CONTINUE
256 l = l1
257 lsv = l
258 lend = m
259 lendsv = lend
260 l1 = m + 1
261 IF( lend.EQ.l )
262 $ GO TO 10
263*
264* Scale submatrix in rows and columns L to LEND
265*
266 anorm = dlanst( 'M', lend-l+1, d( l ), e( l ) )
267 iscale = 0
268 IF( anorm.EQ.zero )
269 $ GO TO 10
270 IF( anorm.GT.ssfmax ) THEN
271 iscale = 1
272 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l+1, 1, d( l ), n,
273 $ info )
274 CALL dlascl( 'G', 0, 0, anorm, ssfmax, lend-l, 1, e( l ), n,
275 $ info )
276 ELSE IF( anorm.LT.ssfmin ) THEN
277 iscale = 2
278 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l+1, 1, d( l ), n,
279 $ info )
280 CALL dlascl( 'G', 0, 0, anorm, ssfmin, lend-l, 1, e( l ), n,
281 $ info )
282 END IF
283*
284* Choose between QL and QR iteration
285*
286 IF( abs( d( lend ) ).LT.abs( d( l ) ) ) THEN
287 lend = lsv
288 l = lendsv
289 END IF
290*
291 IF( lend.GT.l ) THEN
292*
293* QL Iteration
294*
295* Look for small subdiagonal element.
296*
297 40 CONTINUE
298 IF( l.NE.lend ) THEN
299 lendm1 = lend - 1
300 DO 50 m = l, lendm1
301 tst = abs( e( m ) )**2
302 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+
303 $ safmin )GO TO 60
304 50 CONTINUE
305 END IF
306*
307 m = lend
308*
309 60 CONTINUE
310 IF( m.LT.lend )
311 $ e( m ) = zero
312 p = d( l )
313 IF( m.EQ.l )
314 $ GO TO 80
315*
316* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
317* to compute its eigensystem.
318*
319 IF( m.EQ.l+1 ) THEN
320 IF( icompz.GT.0 ) THEN
321 CALL dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )
322 work( l ) = c
323 work( n-1+l ) = s
324 CALL dlasr( 'R', 'V', 'B', n, 2, work( l ),
325 $ work( n-1+l ), z( 1, l ), ldz )
326 ELSE
327 CALL dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )
328 END IF
329 d( l ) = rt1
330 d( l+1 ) = rt2
331 e( l ) = zero
332 l = l + 2
333 IF( l.LE.lend )
334 $ GO TO 40
335 GO TO 140
336 END IF
337*
338 IF( jtot.EQ.nmaxit )
339 $ GO TO 140
340 jtot = jtot + 1
341*
342* Form shift.
343*
344 g = ( d( l+1 )-p ) / ( two*e( l ) )
345 r = dlapy2( g, one )
346 g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )
347*
348 s = one
349 c = one
350 p = zero
351*
352* Inner loop
353*
354 mm1 = m - 1
355 DO 70 i = mm1, l, -1
356 f = s*e( i )
357 b = c*e( i )
358 CALL dlartg( g, f, c, s, r )
359 IF( i.NE.m-1 )
360 $ e( i+1 ) = r
361 g = d( i+1 ) - p
362 r = ( d( i )-g )*s + two*c*b
363 p = s*r
364 d( i+1 ) = g + p
365 g = c*r - b
366*
367* If eigenvectors are desired, then save rotations.
368*
369 IF( icompz.GT.0 ) THEN
370 work( i ) = c
371 work( n-1+i ) = -s
372 END IF
373*
374 70 CONTINUE
375*
376* If eigenvectors are desired, then apply saved rotations.
377*
378 IF( icompz.GT.0 ) THEN
379 mm = m - l + 1
380 CALL dlasr( 'R', 'V', 'B', n, mm, work( l ), work( n-1+l ),
381 $ z( 1, l ), ldz )
382 END IF
383*
384 d( l ) = d( l ) - p
385 e( l ) = g
386 GO TO 40
387*
388* Eigenvalue found.
389*
390 80 CONTINUE
391 d( l ) = p
392*
393 l = l + 1
394 IF( l.LE.lend )
395 $ GO TO 40
396 GO TO 140
397*
398 ELSE
399*
400* QR Iteration
401*
402* Look for small superdiagonal element.
403*
404 90 CONTINUE
405 IF( l.NE.lend ) THEN
406 lendp1 = lend + 1
407 DO 100 m = l, lendp1, -1
408 tst = abs( e( m-1 ) )**2
409 IF( tst.LE.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+
410 $ safmin )GO TO 110
411 100 CONTINUE
412 END IF
413*
414 m = lend
415*
416 110 CONTINUE
417 IF( m.GT.lend )
418 $ e( m-1 ) = zero
419 p = d( l )
420 IF( m.EQ.l )
421 $ GO TO 130
422*
423* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
424* to compute its eigensystem.
425*
426 IF( m.EQ.l-1 ) THEN
427 IF( icompz.GT.0 ) THEN
428 CALL dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )
429 work( m ) = c
430 work( n-1+m ) = s
431 CALL dlasr( 'R', 'V', 'F', n, 2, work( m ),
432 $ work( n-1+m ), z( 1, l-1 ), ldz )
433 ELSE
434 CALL dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )
435 END IF
436 d( l-1 ) = rt1
437 d( l ) = rt2
438 e( l-1 ) = zero
439 l = l - 2
440 IF( l.GE.lend )
441 $ GO TO 90
442 GO TO 140
443 END IF
444*
445 IF( jtot.EQ.nmaxit )
446 $ GO TO 140
447 jtot = jtot + 1
448*
449* Form shift.
450*
451 g = ( d( l-1 )-p ) / ( two*e( l-1 ) )
452 r = dlapy2( g, one )
453 g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )
454*
455 s = one
456 c = one
457 p = zero
458*
459* Inner loop
460*
461 lm1 = l - 1
462 DO 120 i = m, lm1
463 f = s*e( i )
464 b = c*e( i )
465 CALL dlartg( g, f, c, s, r )
466 IF( i.NE.m )
467 $ e( i-1 ) = r
468 g = d( i ) - p
469 r = ( d( i+1 )-g )*s + two*c*b
470 p = s*r
471 d( i ) = g + p
472 g = c*r - b
473*
474* If eigenvectors are desired, then save rotations.
475*
476 IF( icompz.GT.0 ) THEN
477 work( i ) = c
478 work( n-1+i ) = s
479 END IF
480*
481 120 CONTINUE
482*
483* If eigenvectors are desired, then apply saved rotations.
484*
485 IF( icompz.GT.0 ) THEN
486 mm = l - m + 1
487 CALL dlasr( 'R', 'V', 'F', n, mm, work( m ), work( n-1+m ),
488 $ z( 1, m ), ldz )
489 END IF
490*
491 d( l ) = d( l ) - p
492 e( lm1 ) = g
493 GO TO 90
494*
495* Eigenvalue found.
496*
497 130 CONTINUE
498 d( l ) = p
499*
500 l = l - 1
501 IF( l.GE.lend )
502 $ GO TO 90
503 GO TO 140
504*
505 END IF
506*
507* Undo scaling if necessary
508*
509 140 CONTINUE
510 IF( iscale.EQ.1 ) THEN
511 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,
512 $ d( lsv ), n, info )
513 CALL dlascl( 'G', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),
514 $ n, info )
515 ELSE IF( iscale.EQ.2 ) THEN
516 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,
517 $ d( lsv ), n, info )
518 CALL dlascl( 'G', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),
519 $ n, info )
520 END IF
521*
522* Check for no convergence to an eigenvalue after a total
523* of N*MAXIT iterations.
524*
525 IF( jtot.LT.nmaxit )
526 $ GO TO 10
527 DO 150 i = 1, n - 1
528 IF( e( i ).NE.zero )
529 $ info = info + 1
530 150 CONTINUE
531 GO TO 190
532*
533* Order eigenvalues and eigenvectors.
534*
535 160 CONTINUE
536 IF( icompz.EQ.0 ) THEN
537*
538* Use Quick Sort
539*
540 CALL dlasrt( 'I', n, d, info )
541*
542 ELSE
543*
544* Use Selection Sort to minimize swaps of eigenvectors
545*
546 DO 180 ii = 2, n
547 i = ii - 1
548 k = i
549 p = d( i )
550 DO 170 j = ii, n
551 IF( d( j ).LT.p ) THEN
552 k = j
553 p = d( j )
554 END IF
555 170 CONTINUE
556 IF( k.NE.i ) THEN
557 d( k ) = d( i )
558 d( i ) = p
559 CALL dswap( n, z( 1, i ), 1, z( 1, k ), 1 )
560 END IF
561 180 CONTINUE
562 END IF
563*
564 190 CONTINUE
565 RETURN
566*
567* End of DSTEQR
568*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlae2
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
dlaev2
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlanst
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
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlascl
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
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dlasr
subroutine dlasr(side, pivot, direct, m, n, c, s, a, lda)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition dlasr.f:199
dlasrt
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: