LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ssterf.f
Go to the documentation of this file.
00001 *> \brief \b SSTERF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SSTERF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssterf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssterf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssterf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SSTERF( N, D, E, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       REAL               D( * ), E( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> SSTERF 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 REAL 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 REAL 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 SSTERF( 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       REAL               D( * ), E( * )
00099 *     ..
00100 *
00101 *  =====================================================================
00102 *
00103 *     .. Parameters ..
00104       REAL               ZERO, ONE, TWO, THREE
00105       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
00106      $                   THREE = 3.0E0 )
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       REAL               ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
00114      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
00115      $                   SIGMA, SSFMAX, SSFMIN
00116 *     ..
00117 *     .. External Functions ..
00118       REAL               SLAMCH, SLANST, SLAPY2
00119       EXTERNAL           SLAMCH, SLANST, SLAPY2
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           SLAE2, SLASCL, SLASRT, 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( 'SSTERF', -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 = SLAMCH( 'E' )
00146       EPS2 = EPS**2
00147       SAFMIN = SLAMCH( 'S' )
00148       SAFMAX = ONE / SAFMIN
00149       SSFMAX = SQRT( SAFMAX ) / THREE
00150       SSFMIN = SQRT( SAFMIN ) / EPS2
00151 *
00152 *     Compute the eigenvalues of the tridiagonal matrix.
00153 *
00154       NMAXIT = N*MAXIT
00155       SIGMA = ZERO
00156       JTOT = 0
00157 *
00158 *     Determine where the matrix splits and choose QL or QR iteration
00159 *     for each block, according to whether top or bottom diagonal
00160 *     element is smaller.
00161 *
00162       L1 = 1
00163 *
00164    10 CONTINUE
00165       IF( L1.GT.N )
00166      $   GO TO 170
00167       IF( L1.GT.1 )
00168      $   E( L1-1 ) = ZERO
00169       DO 20 M = L1, N - 1
00170          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*
00171      $       SQRT( ABS( D( M+1 ) ) ) )*EPS ) THEN
00172             E( M ) = ZERO
00173             GO TO 30
00174          END IF
00175    20 CONTINUE
00176       M = N
00177 *
00178    30 CONTINUE
00179       L = L1
00180       LSV = L
00181       LEND = M
00182       LENDSV = LEND
00183       L1 = M + 1
00184       IF( LEND.EQ.L )
00185      $   GO TO 10
00186 *
00187 *     Scale submatrix in rows and columns L to LEND
00188 *
00189       ANORM = SLANST( 'M', LEND-L+1, D( L ), E( L ) )
00190       ISCALE = 0
00191       IF( ANORM.EQ.ZERO )
00192      $   GO TO 10      
00193       IF( ANORM.GT.SSFMAX ) THEN
00194          ISCALE = 1
00195          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
00196      $                INFO )
00197          CALL SLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
00198      $                INFO )
00199       ELSE IF( ANORM.LT.SSFMIN ) THEN
00200          ISCALE = 2
00201          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
00202      $                INFO )
00203          CALL SLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
00204      $                INFO )
00205       END IF
00206 *
00207       DO 40 I = L, LEND - 1
00208          E( I ) = E( I )**2
00209    40 CONTINUE
00210 *
00211 *     Choose between QL and QR iteration
00212 *
00213       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
00214          LEND = LSV
00215          L = LENDSV
00216       END IF
00217 *
00218       IF( LEND.GE.L ) THEN
00219 *
00220 *        QL Iteration
00221 *
00222 *        Look for small subdiagonal element.
00223 *
00224    50    CONTINUE
00225          IF( L.NE.LEND ) THEN
00226             DO 60 M = L, LEND - 1
00227                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
00228      $            GO TO 70
00229    60       CONTINUE
00230          END IF
00231          M = LEND
00232 *
00233    70    CONTINUE
00234          IF( M.LT.LEND )
00235      $      E( M ) = ZERO
00236          P = D( L )
00237          IF( M.EQ.L )
00238      $      GO TO 90
00239 *
00240 *        If remaining matrix is 2 by 2, use SLAE2 to compute its
00241 *        eigenvalues.
00242 *
00243          IF( M.EQ.L+1 ) THEN
00244             RTE = SQRT( E( L ) )
00245             CALL SLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
00246             D( L ) = RT1
00247             D( L+1 ) = RT2
00248             E( L ) = ZERO
00249             L = L + 2
00250             IF( L.LE.LEND )
00251      $         GO TO 50
00252             GO TO 150
00253          END IF
00254 *
00255          IF( JTOT.EQ.NMAXIT )
00256      $      GO TO 150
00257          JTOT = JTOT + 1
00258 *
00259 *        Form shift.
00260 *
00261          RTE = SQRT( E( L ) )
00262          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
00263          R = SLAPY2( SIGMA, ONE )
00264          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
00265 *
00266          C = ONE
00267          S = ZERO
00268          GAMMA = D( M ) - SIGMA
00269          P = GAMMA*GAMMA
00270 *
00271 *        Inner loop
00272 *
00273          DO 80 I = M - 1, L, -1
00274             BB = E( I )
00275             R = P + BB
00276             IF( I.NE.M-1 )
00277      $         E( I+1 ) = S*R
00278             OLDC = C
00279             C = P / R
00280             S = BB / R
00281             OLDGAM = GAMMA
00282             ALPHA = D( I )
00283             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
00284             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
00285             IF( C.NE.ZERO ) THEN
00286                P = ( GAMMA*GAMMA ) / C
00287             ELSE
00288                P = OLDC*BB
00289             END IF
00290    80    CONTINUE
00291 *
00292          E( L ) = S*P
00293          D( L ) = SIGMA + GAMMA
00294          GO TO 50
00295 *
00296 *        Eigenvalue found.
00297 *
00298    90    CONTINUE
00299          D( L ) = P
00300 *
00301          L = L + 1
00302          IF( L.LE.LEND )
00303      $      GO TO 50
00304          GO TO 150
00305 *
00306       ELSE
00307 *
00308 *        QR Iteration
00309 *
00310 *        Look for small superdiagonal element.
00311 *
00312   100    CONTINUE
00313          DO 110 M = L, LEND + 1, -1
00314             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
00315      $         GO TO 120
00316   110    CONTINUE
00317          M = LEND
00318 *
00319   120    CONTINUE
00320          IF( M.GT.LEND )
00321      $      E( M-1 ) = ZERO
00322          P = D( L )
00323          IF( M.EQ.L )
00324      $      GO TO 140
00325 *
00326 *        If remaining matrix is 2 by 2, use SLAE2 to compute its
00327 *        eigenvalues.
00328 *
00329          IF( M.EQ.L-1 ) THEN
00330             RTE = SQRT( E( L-1 ) )
00331             CALL SLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
00332             D( L ) = RT1
00333             D( L-1 ) = RT2
00334             E( L-1 ) = ZERO
00335             L = L - 2
00336             IF( L.GE.LEND )
00337      $         GO TO 100
00338             GO TO 150
00339          END IF
00340 *
00341          IF( JTOT.EQ.NMAXIT )
00342      $      GO TO 150
00343          JTOT = JTOT + 1
00344 *
00345 *        Form shift.
00346 *
00347          RTE = SQRT( E( L-1 ) )
00348          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
00349          R = SLAPY2( SIGMA, ONE )
00350          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
00351 *
00352          C = ONE
00353          S = ZERO
00354          GAMMA = D( M ) - SIGMA
00355          P = GAMMA*GAMMA
00356 *
00357 *        Inner loop
00358 *
00359          DO 130 I = M, L - 1
00360             BB = E( I )
00361             R = P + BB
00362             IF( I.NE.M )
00363      $         E( I-1 ) = S*R
00364             OLDC = C
00365             C = P / R
00366             S = BB / R
00367             OLDGAM = GAMMA
00368             ALPHA = D( I+1 )
00369             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
00370             D( I ) = OLDGAM + ( ALPHA-GAMMA )
00371             IF( C.NE.ZERO ) THEN
00372                P = ( GAMMA*GAMMA ) / C
00373             ELSE
00374                P = OLDC*BB
00375             END IF
00376   130    CONTINUE
00377 *
00378          E( L-1 ) = S*P
00379          D( L ) = SIGMA + GAMMA
00380          GO TO 100
00381 *
00382 *        Eigenvalue found.
00383 *
00384   140    CONTINUE
00385          D( L ) = P
00386 *
00387          L = L - 1
00388          IF( L.GE.LEND )
00389      $      GO TO 100
00390          GO TO 150
00391 *
00392       END IF
00393 *
00394 *     Undo scaling if necessary
00395 *
00396   150 CONTINUE
00397       IF( ISCALE.EQ.1 )
00398      $   CALL SLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
00399      $                D( LSV ), N, INFO )
00400       IF( ISCALE.EQ.2 )
00401      $   CALL SLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
00402      $                D( LSV ), N, INFO )
00403 *
00404 *     Check for no convergence to an eigenvalue after a total
00405 *     of N*MAXIT iterations.
00406 *
00407       IF( JTOT.LT.NMAXIT )
00408      $   GO TO 10
00409       DO 160 I = 1, N - 1
00410          IF( E( I ).NE.ZERO )
00411      $      INFO = INFO + 1
00412   160 CONTINUE
00413       GO TO 180
00414 *
00415 *     Sort eigenvalues in increasing order.
00416 *
00417   170 CONTINUE
00418       CALL SLASRT( 'I', N, D, INFO )
00419 *
00420   180 CONTINUE
00421       RETURN
00422 *
00423 *     End of SSTERF
00424 *
00425       END
 All Files Functions