427 lines
10 KiB
Fortran
427 lines
10 KiB
Fortran
*> \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 November 2011
|
|
*
|
|
*> \ingroup auxOTHERcomputational
|
|
*
|
|
* =====================================================================
|
|
SUBROUTINE DSTERF( N, D, E, INFO )
|
|
*
|
|
* -- LAPACK computational routine (version 3.4.0) --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
* November 2011
|
|
*
|
|
* .. 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
|