LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dsterf.f
Go to the documentation of this file.
00001 *> \brief \b DSTERF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DSTERF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DSTERF( N, D, E, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       DOUBLE PRECISION   D( * ), E( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
00037 *> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
00038 *> \endverbatim
00039 *
00040 *  Arguments:
00041 *  ==========
00042 *
00043 *> \param[in] N
00044 *> \verbatim
00045 *>          N is INTEGER
00046 *>          The order of the matrix.  N >= 0.
00047 *> \endverbatim
00048 *>
00049 *> \param[in,out] D
00050 *> \verbatim
00051 *>          D is DOUBLE PRECISION array, dimension (N)
00052 *>          On entry, the n diagonal elements of the tridiagonal matrix.
00053 *>          On exit, if INFO = 0, the eigenvalues in ascending order.
00054 *> \endverbatim
00055 *>
00056 *> \param[in,out] E
00057 *> \verbatim
00058 *>          E is DOUBLE PRECISION array, dimension (N-1)
00059 *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
00060 *>          matrix.
00061 *>          On exit, E has been destroyed.
00062 *> \endverbatim
00063 *>
00064 *> \param[out] INFO
00065 *> \verbatim
00066 *>          INFO is INTEGER
00067 *>          = 0:  successful exit
00068 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00069 *>          > 0:  the algorithm failed to find all of the eigenvalues in
00070 *>                a total of 30*N iterations; if INFO = i, then i
00071 *>                elements of E have not converged to zero.
00072 *> \endverbatim
00073 *
00074 *  Authors:
00075 *  ========
00076 *
00077 *> \author Univ. of Tennessee 
00078 *> \author Univ. of California Berkeley 
00079 *> \author Univ. of Colorado Denver 
00080 *> \author NAG Ltd. 
00081 *
00082 *> \date November 2011
00083 *
00084 *> \ingroup auxOTHERcomputational
00085 *
00086 *  =====================================================================
00087       SUBROUTINE DSTERF( N, D, E, INFO )
00088 *
00089 *  -- LAPACK computational routine (version 3.4.0) --
00090 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00091 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00092 *     November 2011
00093 *
00094 *     .. Scalar Arguments ..
00095       INTEGER            INFO, N
00096 *     ..
00097 *     .. Array Arguments ..
00098       DOUBLE PRECISION   D( * ), E( * )
00099 *     ..
00100 *
00101 *  =====================================================================
00102 *
00103 *     .. Parameters ..
00104       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
00105       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
00106      $                   THREE = 3.0D0 )
00107       INTEGER            MAXIT
00108       PARAMETER          ( MAXIT = 30 )
00109 *     ..
00110 *     .. Local Scalars ..
00111       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
00112      $                   NMAXIT
00113       DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
00114      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
00115      $                   SIGMA, SSFMAX, SSFMIN, RMAX
00116 *     ..
00117 *     .. External Functions ..
00118       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
00119       EXTERNAL           DLAMCH, DLANST, DLAPY2
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           DLAE2, DLASCL, DLASRT, XERBLA
00123 *     ..
00124 *     .. Intrinsic Functions ..
00125       INTRINSIC          ABS, SIGN, SQRT
00126 *     ..
00127 *     .. Executable Statements ..
00128 *
00129 *     Test the input parameters.
00130 *
00131       INFO = 0
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( N.LT.0 ) THEN
00136          INFO = -1
00137          CALL XERBLA( 'DSTERF', -INFO )
00138          RETURN
00139       END IF
00140       IF( N.LE.1 )
00141      $   RETURN
00142 *
00143 *     Determine the unit roundoff for this environment.
00144 *
00145       EPS = DLAMCH( 'E' )
00146       EPS2 = EPS**2
00147       SAFMIN = DLAMCH( 'S' )
00148       SAFMAX = ONE / SAFMIN
00149       SSFMAX = SQRT( SAFMAX ) / THREE
00150       SSFMIN = SQRT( SAFMIN ) / EPS2
00151       RMAX = DLAMCH( 'O' )
00152 *
00153 *     Compute the eigenvalues of the tridiagonal matrix.
00154 *
00155       NMAXIT = N*MAXIT
00156       SIGMA = ZERO
00157       JTOT = 0
00158 *
00159 *     Determine where the matrix splits and choose QL or QR iteration
00160 *     for each block, according to whether top or bottom diagonal
00161 *     element is smaller.
00162 *
00163       L1 = 1
00164 *
00165    10 CONTINUE
00166       IF( L1.GT.N )
00167      $   GO TO 170
00168       IF( L1.GT.1 )
00169      $   E( L1-1 ) = ZERO
00170       DO 20 M = L1, N - 1
00171          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
00172      $       1 ) ) ) )*EPS ) THEN
00173             E( M ) = ZERO
00174             GO TO 30
00175          END IF
00176    20 CONTINUE
00177       M = N
00178 *
00179    30 CONTINUE
00180       L = L1
00181       LSV = L
00182       LEND = M
00183       LENDSV = LEND
00184       L1 = M + 1
00185       IF( LEND.EQ.L )
00186      $   GO TO 10
00187 *
00188 *     Scale submatrix in rows and columns L to LEND
00189 *
00190       ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
00191       ISCALE = 0
00192       IF( ANORM.EQ.ZERO )
00193      $   GO TO 10      
00194       IF( (ANORM.GT.SSFMAX) ) THEN
00195          ISCALE = 1
00196          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
00197      $                INFO )
00198          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
00199      $                INFO )
00200       ELSE IF( ANORM.LT.SSFMIN ) THEN
00201          ISCALE = 2
00202          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
00203      $                INFO )
00204          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
00205      $                INFO )
00206       END IF
00207 *
00208       DO 40 I = L, LEND - 1
00209          E( I ) = E( I )**2
00210    40 CONTINUE
00211 *
00212 *     Choose between QL and QR iteration
00213 *
00214       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
00215          LEND = LSV
00216          L = LENDSV
00217       END IF
00218 *
00219       IF( LEND.GE.L ) THEN
00220 *
00221 *        QL Iteration
00222 *
00223 *        Look for small subdiagonal element.
00224 *
00225    50    CONTINUE
00226          IF( L.NE.LEND ) THEN
00227             DO 60 M = L, LEND - 1
00228                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
00229      $            GO TO 70
00230    60       CONTINUE
00231          END IF
00232          M = LEND
00233 *
00234    70    CONTINUE
00235          IF( M.LT.LEND )
00236      $      E( M ) = ZERO
00237          P = D( L )
00238          IF( M.EQ.L )
00239      $      GO TO 90
00240 *
00241 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
00242 *        eigenvalues.
00243 *
00244          IF( M.EQ.L+1 ) THEN
00245             RTE = SQRT( E( L ) )
00246             CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
00247             D( L ) = RT1
00248             D( L+1 ) = RT2
00249             E( L ) = ZERO
00250             L = L + 2
00251             IF( L.LE.LEND )
00252      $         GO TO 50
00253             GO TO 150
00254          END IF
00255 *
00256          IF( JTOT.EQ.NMAXIT )
00257      $      GO TO 150
00258          JTOT = JTOT + 1
00259 *
00260 *        Form shift.
00261 *
00262          RTE = SQRT( E( L ) )
00263          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
00264          R = DLAPY2( SIGMA, ONE )
00265          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
00266 *
00267          C = ONE
00268          S = ZERO
00269          GAMMA = D( M ) - SIGMA
00270          P = GAMMA*GAMMA
00271 *
00272 *        Inner loop
00273 *
00274          DO 80 I = M - 1, L, -1
00275             BB = E( I )
00276             R = P + BB
00277             IF( I.NE.M-1 )
00278      $         E( I+1 ) = S*R
00279             OLDC = C
00280             C = P / R
00281             S = BB / R
00282             OLDGAM = GAMMA
00283             ALPHA = D( I )
00284             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
00285             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
00286             IF( C.NE.ZERO ) THEN
00287                P = ( GAMMA*GAMMA ) / C
00288             ELSE
00289                P = OLDC*BB
00290             END IF
00291    80    CONTINUE
00292 *
00293          E( L ) = S*P
00294          D( L ) = SIGMA + GAMMA
00295          GO TO 50
00296 *
00297 *        Eigenvalue found.
00298 *
00299    90    CONTINUE
00300          D( L ) = P
00301 *
00302          L = L + 1
00303          IF( L.LE.LEND )
00304      $      GO TO 50
00305          GO TO 150
00306 *
00307       ELSE
00308 *
00309 *        QR Iteration
00310 *
00311 *        Look for small superdiagonal element.
00312 *
00313   100    CONTINUE
00314          DO 110 M = L, LEND + 1, -1
00315             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
00316      $         GO TO 120
00317   110    CONTINUE
00318          M = LEND
00319 *
00320   120    CONTINUE
00321          IF( M.GT.LEND )
00322      $      E( M-1 ) = ZERO
00323          P = D( L )
00324          IF( M.EQ.L )
00325      $      GO TO 140
00326 *
00327 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
00328 *        eigenvalues.
00329 *
00330          IF( M.EQ.L-1 ) THEN
00331             RTE = SQRT( E( L-1 ) )
00332             CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
00333             D( L ) = RT1
00334             D( L-1 ) = RT2
00335             E( L-1 ) = ZERO
00336             L = L - 2
00337             IF( L.GE.LEND )
00338      $         GO TO 100
00339             GO TO 150
00340          END IF
00341 *
00342          IF( JTOT.EQ.NMAXIT )
00343      $      GO TO 150
00344          JTOT = JTOT + 1
00345 *
00346 *        Form shift.
00347 *
00348          RTE = SQRT( E( L-1 ) )
00349          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
00350          R = DLAPY2( SIGMA, ONE )
00351          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
00352 *
00353          C = ONE
00354          S = ZERO
00355          GAMMA = D( M ) - SIGMA
00356          P = GAMMA*GAMMA
00357 *
00358 *        Inner loop
00359 *
00360          DO 130 I = M, L - 1
00361             BB = E( I )
00362             R = P + BB
00363             IF( I.NE.M )
00364      $         E( I-1 ) = S*R
00365             OLDC = C
00366             C = P / R
00367             S = BB / R
00368             OLDGAM = GAMMA
00369             ALPHA = D( I+1 )
00370             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
00371             D( I ) = OLDGAM + ( ALPHA-GAMMA )
00372             IF( C.NE.ZERO ) THEN
00373                P = ( GAMMA*GAMMA ) / C
00374             ELSE
00375                P = OLDC*BB
00376             END IF
00377   130    CONTINUE
00378 *
00379          E( L-1 ) = S*P
00380          D( L ) = SIGMA + GAMMA
00381          GO TO 100
00382 *
00383 *        Eigenvalue found.
00384 *
00385   140    CONTINUE
00386          D( L ) = P
00387 *
00388          L = L - 1
00389          IF( L.GE.LEND )
00390      $      GO TO 100
00391          GO TO 150
00392 *
00393       END IF
00394 *
00395 *     Undo scaling if necessary
00396 *
00397   150 CONTINUE
00398       IF( ISCALE.EQ.1 )
00399      $   CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
00400      $                D( LSV ), N, INFO )
00401       IF( ISCALE.EQ.2 )
00402      $   CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
00403      $                D( LSV ), N, INFO )
00404 *
00405 *     Check for no convergence to an eigenvalue after a total
00406 *     of N*MAXIT iterations.
00407 *
00408       IF( JTOT.LT.NMAXIT )
00409      $   GO TO 10
00410       DO 160 I = 1, N - 1
00411          IF( E( I ).NE.ZERO )
00412      $      INFO = INFO + 1
00413   160 CONTINUE
00414       GO TO 180
00415 *
00416 *     Sort eigenvalues in increasing order.
00417 *
00418   170 CONTINUE
00419       CALL DLASRT( 'I', N, D, INFO )
00420 *
00421   180 CONTINUE
00422       RETURN
00423 *
00424 *     End of DSTERF
00425 *
00426       END
 All Files Functions