LAPACK 3.3.0

dlarrk.f

Go to the documentation of this file.
00001       SUBROUTINE DLARRK( N, IW, GL, GU,
00002      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
00003       IMPLICIT NONE
00004 *
00005 *  -- LAPACK auxiliary routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER   INFO, IW, N
00012       DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   D( * ), E2( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLARRK computes one eigenvalue of a symmetric tridiagonal
00022 *  matrix T to suitable accuracy. This is an auxiliary code to be
00023 *  called from DSTEMR.
00024 *
00025 *  To avoid overflow, the matrix must be scaled so that its
00026 *  largest element is no greater than overflow**(1/2) *
00027 *  underflow**(1/4) in absolute value, and for greatest
00028 *  accuracy, it should not be much smaller than that.
00029 *
00030 *  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00031 *  Matrix", Report CS41, Computer Science Dept., Stanford
00032 *  University, July 21, 1966.
00033 *
00034 *  Arguments
00035 *  =========
00036 *
00037 *  N       (input) INTEGER
00038 *          The order of the tridiagonal matrix T.  N >= 0.
00039 *
00040 *  IW      (input) INTEGER
00041 *          The index of the eigenvalues to be returned.
00042 *
00043 *  GL      (input) DOUBLE PRECISION
00044 *  GU      (input) DOUBLE PRECISION
00045 *          An upper and a lower bound on the eigenvalue.
00046 *
00047 *  D       (input) DOUBLE PRECISION array, dimension (N)
00048 *          The n diagonal elements of the tridiagonal matrix T.
00049 *
00050 *  E2      (input) DOUBLE PRECISION array, dimension (N-1)
00051 *          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
00052 *
00053 *  PIVMIN  (input) DOUBLE PRECISION
00054 *          The minimum pivot allowed in the Sturm sequence for T.
00055 *
00056 *  RELTOL  (input) DOUBLE PRECISION
00057 *          The minimum relative width of an interval.  When an interval
00058 *          is narrower than RELTOL times the larger (in
00059 *          magnitude) endpoint, then it is considered to be
00060 *          sufficiently small, i.e., converged.  Note: this should
00061 *          always be at least radix*machine epsilon.
00062 *
00063 *  W       (output) DOUBLE PRECISION
00064 *
00065 *  WERR    (output) DOUBLE PRECISION
00066 *          The error bound on the corresponding eigenvalue approximation
00067 *          in W.
00068 *
00069 *  INFO    (output) INTEGER
00070 *          = 0:       Eigenvalue converged
00071 *          = -1:      Eigenvalue did NOT converge
00072 *
00073 *  Internal Parameters
00074 *  ===================
00075 *
00076 *  FUDGE   DOUBLE PRECISION, default = 2
00077 *          A "fudge factor" to widen the Gershgorin intervals.
00078 *
00079 *  =====================================================================
00080 *
00081 *     .. Parameters ..
00082       DOUBLE PRECISION   FUDGE, HALF, TWO, ZERO
00083       PARAMETER          ( HALF = 0.5D0, TWO = 2.0D0,
00084      $                     FUDGE = TWO, ZERO = 0.0D0 )
00085 *     ..
00086 *     .. Local Scalars ..
00087       INTEGER   I, IT, ITMAX, NEGCNT
00088       DOUBLE PRECISION   ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
00089      $                   TMP2, TNORM
00090 *     ..
00091 *     .. External Functions ..
00092       DOUBLE PRECISION   DLAMCH
00093       EXTERNAL   DLAMCH
00094 *     ..
00095 *     .. Intrinsic Functions ..
00096       INTRINSIC          ABS, INT, LOG, MAX
00097 *     ..
00098 *     .. Executable Statements ..
00099 *
00100 *     Get machine constants
00101       EPS = DLAMCH( 'P' )
00102 
00103       TNORM = MAX( ABS( GL ), ABS( GU ) )
00104       RTOLI = RELTOL
00105       ATOLI = FUDGE*TWO*PIVMIN
00106 
00107       ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
00108      $           LOG( TWO ) ) + 2
00109 
00110       INFO = -1
00111 
00112       LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
00113       RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
00114       IT = 0
00115 
00116  10   CONTINUE
00117 *
00118 *     Check if interval converged or maximum number of iterations reached
00119 *
00120       TMP1 = ABS( RIGHT - LEFT )
00121       TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
00122       IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
00123          INFO = 0
00124          GOTO 30
00125       ENDIF
00126       IF(IT.GT.ITMAX)
00127      $   GOTO 30
00128 
00129 *
00130 *     Count number of negative pivots for mid-point
00131 *
00132       IT = IT + 1
00133       MID = HALF * (LEFT + RIGHT)
00134       NEGCNT = 0
00135       TMP1 = D( 1 ) - MID
00136       IF( ABS( TMP1 ).LT.PIVMIN )
00137      $   TMP1 = -PIVMIN
00138       IF( TMP1.LE.ZERO )
00139      $   NEGCNT = NEGCNT + 1
00140 *
00141       DO 20 I = 2, N
00142          TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
00143          IF( ABS( TMP1 ).LT.PIVMIN )
00144      $      TMP1 = -PIVMIN
00145          IF( TMP1.LE.ZERO )
00146      $      NEGCNT = NEGCNT + 1
00147  20   CONTINUE
00148 
00149       IF(NEGCNT.GE.IW) THEN
00150          RIGHT = MID
00151       ELSE
00152          LEFT = MID
00153       ENDIF
00154       GOTO 10
00155 
00156  30   CONTINUE
00157 *
00158 *     Converged or maximum number of iterations reached
00159 *
00160       W = HALF * (LEFT + RIGHT)
00161       WERR = HALF * ABS( RIGHT - LEFT )
00162 
00163       RETURN
00164 *
00165 *     End of DLARRK
00166 *
00167       END
 All Files Functions