#include "ESMF_LapackBlas.inc" !> \brief \b DSTERF ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DSTERF + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DSTERF( N, D, E, INFO ) ! ! .. Scalar Arguments .. ! INTEGER INFO, N ! .. ! .. Array Arguments .. ! DOUBLE PRECISION D( * ), E( * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix !> using the Pal-Walker-Kahan variant of the QL or QR algorithm. !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] N !> \verbatim !> N is INTEGER !> The order of the matrix. N >= 0. !> \endverbatim !> !> \param[in,out] D !> \verbatim !> D is DOUBLE PRECISION array, dimension (N) !> On entry, the n diagonal elements of the tridiagonal matrix. !> On exit, if INFO = 0, the eigenvalues in ascending order. !> \endverbatim !> !> \param[in,out] E !> \verbatim !> 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. !> \endverbatim !> !> \param[out] INFO !> \verbatim !> INFO is INTEGER !> = 0: successful exit !> < 0: if INFO = -i, the i-th argument had an illegal value !> > 0: the algorithm failed to find all of the eigenvalues in !> a total of 30*N iterations; if INFO = i, then i !> elements of E have not converged to zero. !> \endverbatim ! ! Authors: ! ======== ! !> \author Univ. of Tennessee !> \author Univ. of California Berkeley !> \author Univ. of Colorado Denver !> \author NAG Ltd. ! !> \date December 2016 ! !> \ingroup auxOTHERcomputational ! ! ===================================================================== SUBROUTINE DSTERF( N, D, E, INFO ) ! ! -- LAPACK computational routine (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. Scalar Arguments .. INTEGER INFO, N ! .. ! .. Array Arguments .. DOUBLE PRECISION D( * ), E( * ) ! .. ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, THREE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, & THREE = 3.0D0 ) INTEGER MAXIT PARAMETER ( MAXIT = 30 ) ! .. ! .. Local Scalars .. INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M, & NMAXIT DOUBLE PRECISION ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC, & OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN, & SIGMA, SSFMAX, SSFMIN, RMAX ! .. ! .. External Functions .. DOUBLE PRECISION DLAMCH, DLANST, DLAPY2 EXTERNAL DLAMCH, DLANST, DLAPY2 ! .. ! .. External Subroutines .. EXTERNAL DLAE2, DLASCL, DLASRT, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, SIGN, SQRT ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 ! ! Quick return if possible ! IF( N.LT.0 ) THEN INFO = -1 CALL XERBLA( 'DSTERF', -INFO ) RETURN END IF IF( N.LE.1 ) & RETURN ! ! Determine the unit roundoff for this environment. ! EPS = DLAMCH( 'E' ) EPS2 = EPS**2 SAFMIN = DLAMCH( 'S' ) SAFMAX = ONE / SAFMIN SSFMAX = SQRT( SAFMAX ) / THREE SSFMIN = SQRT( SAFMIN ) / EPS2 RMAX = DLAMCH( 'O' ) ! ! Compute the eigenvalues of the tridiagonal matrix. ! NMAXIT = N*MAXIT SIGMA = ZERO JTOT = 0 ! ! Determine where the matrix splits and choose QL or QR iteration ! for each block, according to whether top or bottom diagonal ! element is smaller. ! L1 = 1 ! 10 CONTINUE IF( L1.GT.N ) & GO TO 170 IF( L1.GT.1 ) & E( L1-1 ) = ZERO DO 20 M = L1, N - 1 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+ & 1 ) ) ) )*EPS ) THEN E( M ) = ZERO GO TO 30 END IF 20 CONTINUE M = N ! 30 CONTINUE L = L1 LSV = L LEND = M LENDSV = LEND L1 = M + 1 IF( LEND.EQ.L ) & GO TO 10 ! ! Scale submatrix in rows and columns L to LEND ! ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) ISCALE = 0 IF( ANORM.EQ.ZERO ) & GO TO 10 IF( (ANORM.GT.SSFMAX) ) THEN ISCALE = 1 CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N, & INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N, & INFO ) ELSE IF( ANORM.LT.SSFMIN ) THEN ISCALE = 2 CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N, & INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N, & INFO ) END IF ! DO 40 I = L, LEND - 1 E( I ) = E( I )**2 40 CONTINUE ! ! Choose between QL and QR iteration ! IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN LEND = LSV L = LENDSV END IF ! IF( LEND.GE.L ) THEN ! ! QL Iteration ! ! Look for small subdiagonal element. ! 50 CONTINUE IF( L.NE.LEND ) THEN DO 60 M = L, LEND - 1 IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) ) & GO TO 70 60 CONTINUE END IF M = LEND ! 70 CONTINUE IF( M.LT.LEND ) & E( M ) = ZERO P = D( L ) IF( M.EQ.L ) & GO TO 90 ! ! If remaining matrix is 2 by 2, use DLAE2 to compute its ! eigenvalues. ! IF( M.EQ.L+1 ) THEN RTE = SQRT( E( L ) ) CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) D( L ) = RT1 D( L+1 ) = RT2 E( L ) = ZERO L = L + 2 IF( L.LE.LEND ) & GO TO 50 GO TO 150 END IF ! IF( JTOT.EQ.NMAXIT ) & GO TO 150 JTOT = JTOT + 1 ! ! Form shift. ! RTE = SQRT( E( L ) ) SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) R = DLAPY2( SIGMA, ONE ) SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) ! C = ONE S = ZERO GAMMA = D( M ) - SIGMA P = GAMMA*GAMMA ! ! Inner loop ! DO 80 I = M - 1, L, -1 BB = E( I ) R = P + BB IF( I.NE.M-1 ) & E( I+1 ) = S*R OLDC = C C = P / R S = BB / R OLDGAM = GAMMA ALPHA = D( I ) GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) IF( C.NE.ZERO ) THEN P = ( GAMMA*GAMMA ) / C ELSE P = OLDC*BB END IF 80 CONTINUE ! E( L ) = S*P D( L ) = SIGMA + GAMMA GO TO 50 ! ! Eigenvalue found. ! 90 CONTINUE D( L ) = P ! L = L + 1 IF( L.LE.LEND ) & GO TO 50 GO TO 150 ! ELSE ! ! QR Iteration ! ! Look for small superdiagonal element. ! 100 CONTINUE DO 110 M = L, LEND + 1, -1 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) ) & GO TO 120 110 CONTINUE M = LEND ! 120 CONTINUE IF( M.GT.LEND ) & E( M-1 ) = ZERO P = D( L ) IF( M.EQ.L ) & GO TO 140 ! ! If remaining matrix is 2 by 2, use DLAE2 to compute its ! eigenvalues. ! IF( M.EQ.L-1 ) THEN RTE = SQRT( E( L-1 ) ) CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 ) D( L ) = RT1 D( L-1 ) = RT2 E( L-1 ) = ZERO L = L - 2 IF( L.GE.LEND ) & GO TO 100 GO TO 150 END IF ! IF( JTOT.EQ.NMAXIT ) & GO TO 150 JTOT = JTOT + 1 ! ! Form shift. ! RTE = SQRT( E( L-1 ) ) SIGMA = ( D( L-1 )-P ) / ( TWO*RTE ) R = DLAPY2( SIGMA, ONE ) SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) ) ! C = ONE S = ZERO GAMMA = D( M ) - SIGMA P = GAMMA*GAMMA ! ! Inner loop ! DO 130 I = M, L - 1 BB = E( I ) R = P + BB IF( I.NE.M ) & E( I-1 ) = S*R OLDC = C C = P / R S = BB / R OLDGAM = GAMMA ALPHA = D( I+1 ) GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM D( I ) = OLDGAM + ( ALPHA-GAMMA ) IF( C.NE.ZERO ) THEN P = ( GAMMA*GAMMA ) / C ELSE P = OLDC*BB END IF 130 CONTINUE ! E( L-1 ) = S*P D( L ) = SIGMA + GAMMA GO TO 100 ! ! Eigenvalue found. ! 140 CONTINUE D( L ) = P ! L = L - 1 IF( L.GE.LEND ) & GO TO 100 GO TO 150 ! END IF ! ! Undo scaling if necessary ! 150 CONTINUE IF( ISCALE.EQ.1 ) & CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1, & D( LSV ), N, INFO ) IF( ISCALE.EQ.2 ) & CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1, & D( LSV ), N, INFO ) ! ! Check for no convergence to an eigenvalue after a total ! of N*MAXIT iterations. ! IF( JTOT.LT.NMAXIT ) & GO TO 10 DO 160 I = 1, N - 1 IF( E( I ).NE.ZERO ) & INFO = INFO + 1 160 CONTINUE GO TO 180 ! ! Sort eigenvalues in increasing order. ! 170 CONTINUE CALL DLASRT( 'I', N, D, INFO ) ! 180 CONTINUE RETURN ! ! End of DSTERF ! END