| *> \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 |
| print *, "l1 = ", l1 |
| IF( L1.GT.N ) THEN |
| print *, "going to 170" |
| GO TO 170 |
| end if |
| 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 |
| print *, "30, d" |
| print *, d(1:n) |
| 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 |
| print *, "ql, d" |
| print *, d(1:n) |
| * |
| * 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 |
| * |
| print *, "inner loop d before" |
| print *, d(1:n) |
| |
| DO 80 I = M - 1, L, -1 |
| print *, "inner loop" |
| print *, "ei", e(i) |
| BB = E( I ) |
| R = P + BB |
| print *, "bb,p,r" |
| print *, bb,p,r |
| IF( I.NE.M-1 ) THEN |
| print *, s,r |
| E( I+1 ) = S*R |
| end if |
| OLDC = C |
| C = P / R |
| S = BB / R |
| OLDGAM = GAMMA |
| print *, "di", d(i) |
| ALPHA = D( I ) |
| GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM |
| print *,"og, a, ga", OLDGAM, ALPHA, GAMMA |
| D( I+1 ) = OLDGAM + ( ALPHA-GAMMA ) |
| IF( C.NE.ZERO ) THEN |
| P = ( GAMMA*GAMMA ) / C |
| ELSE |
| P = OLDC*BB |
| END IF |
| print *, "p, gamma = ", p,GAMMA |
| 80 CONTINUE |
| * |
| E( L ) = S*P |
| D( L ) = SIGMA + GAMMA |
| |
| print *, "inner loop d after" |
| print *, d(1:n) |
| 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 |
