Annotation of rpl/lapack/lapack/dlarrk.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b DLARRK
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLARRK + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrk.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrk.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrk.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLARRK( N, IW, GL, GU,
! 22: * D, E2, PIVMIN, RELTOL, W, WERR, INFO)
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, IW, N
! 26: * DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION D( * ), E2( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DLARRK computes one eigenvalue of a symmetric tridiagonal
! 39: *> matrix T to suitable accuracy. This is an auxiliary code to be
! 40: *> called from DSTEMR.
! 41: *>
! 42: *> To avoid overflow, the matrix must be scaled so that its
! 43: *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
! 44: *> accuracy, it should not be much smaller than that.
! 45: *>
! 46: *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
! 47: *> Matrix", Report CS41, Computer Science Dept., Stanford
! 48: *> University, July 21, 1966.
! 49: *> \endverbatim
! 50: *
! 51: * Arguments:
! 52: * ==========
! 53: *
! 54: *> \param[in] N
! 55: *> \verbatim
! 56: *> N is INTEGER
! 57: *> The order of the tridiagonal matrix T. N >= 0.
! 58: *> \endverbatim
! 59: *>
! 60: *> \param[in] IW
! 61: *> \verbatim
! 62: *> IW is INTEGER
! 63: *> The index of the eigenvalues to be returned.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] GL
! 67: *> \verbatim
! 68: *> GL is DOUBLE PRECISION
! 69: *> \endverbatim
! 70: *>
! 71: *> \param[in] GU
! 72: *> \verbatim
! 73: *> GU is DOUBLE PRECISION
! 74: *> An upper and a lower bound on the eigenvalue.
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in] D
! 78: *> \verbatim
! 79: *> D is DOUBLE PRECISION array, dimension (N)
! 80: *> The n diagonal elements of the tridiagonal matrix T.
! 81: *> \endverbatim
! 82: *>
! 83: *> \param[in] E2
! 84: *> \verbatim
! 85: *> E2 is DOUBLE PRECISION array, dimension (N-1)
! 86: *> The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] PIVMIN
! 90: *> \verbatim
! 91: *> PIVMIN is DOUBLE PRECISION
! 92: *> The minimum pivot allowed in the Sturm sequence for T.
! 93: *> \endverbatim
! 94: *>
! 95: *> \param[in] RELTOL
! 96: *> \verbatim
! 97: *> RELTOL is DOUBLE PRECISION
! 98: *> The minimum relative width of an interval. When an interval
! 99: *> is narrower than RELTOL times the larger (in
! 100: *> magnitude) endpoint, then it is considered to be
! 101: *> sufficiently small, i.e., converged. Note: this should
! 102: *> always be at least radix*machine epsilon.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[out] W
! 106: *> \verbatim
! 107: *> W is DOUBLE PRECISION
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[out] WERR
! 111: *> \verbatim
! 112: *> WERR is DOUBLE PRECISION
! 113: *> The error bound on the corresponding eigenvalue approximation
! 114: *> in W.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[out] INFO
! 118: *> \verbatim
! 119: *> INFO is INTEGER
! 120: *> = 0: Eigenvalue converged
! 121: *> = -1: Eigenvalue did NOT converge
! 122: *> \endverbatim
! 123: *
! 124: *> \par Internal Parameters:
! 125: * =========================
! 126: *>
! 127: *> \verbatim
! 128: *> FUDGE DOUBLE PRECISION, default = 2
! 129: *> A "fudge factor" to widen the Gershgorin intervals.
! 130: *> \endverbatim
! 131: *
! 132: * Authors:
! 133: * ========
! 134: *
! 135: *> \author Univ. of Tennessee
! 136: *> \author Univ. of California Berkeley
! 137: *> \author Univ. of Colorado Denver
! 138: *> \author NAG Ltd.
! 139: *
! 140: *> \date November 2011
! 141: *
! 142: *> \ingroup auxOTHERauxiliary
! 143: *
! 144: * =====================================================================
1.1 bertrand 145: SUBROUTINE DLARRK( N, IW, GL, GU,
146: $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
147: *
1.8 ! bertrand 148: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 149: * -- LAPACK is a software package provided by Univ. of Tennessee, --
150: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 151: * November 2011
1.1 bertrand 152: *
153: * .. Scalar Arguments ..
154: INTEGER INFO, IW, N
155: DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
156: * ..
157: * .. Array Arguments ..
158: DOUBLE PRECISION D( * ), E2( * )
159: * ..
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
165: PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
166: $ FUDGE = TWO, ZERO = 0.0D0 )
167: * ..
168: * .. Local Scalars ..
169: INTEGER I, IT, ITMAX, NEGCNT
170: DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
171: $ TMP2, TNORM
172: * ..
173: * .. External Functions ..
174: DOUBLE PRECISION DLAMCH
175: EXTERNAL DLAMCH
176: * ..
177: * .. Intrinsic Functions ..
178: INTRINSIC ABS, INT, LOG, MAX
179: * ..
180: * .. Executable Statements ..
181: *
182: * Get machine constants
183: EPS = DLAMCH( 'P' )
184:
185: TNORM = MAX( ABS( GL ), ABS( GU ) )
186: RTOLI = RELTOL
187: ATOLI = FUDGE*TWO*PIVMIN
188:
189: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
190: $ LOG( TWO ) ) + 2
191:
192: INFO = -1
193:
194: LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
195: RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
196: IT = 0
197:
198: 10 CONTINUE
199: *
200: * Check if interval converged or maximum number of iterations reached
201: *
202: TMP1 = ABS( RIGHT - LEFT )
203: TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
204: IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
205: INFO = 0
206: GOTO 30
207: ENDIF
208: IF(IT.GT.ITMAX)
209: $ GOTO 30
210:
211: *
212: * Count number of negative pivots for mid-point
213: *
214: IT = IT + 1
215: MID = HALF * (LEFT + RIGHT)
216: NEGCNT = 0
217: TMP1 = D( 1 ) - MID
218: IF( ABS( TMP1 ).LT.PIVMIN )
219: $ TMP1 = -PIVMIN
220: IF( TMP1.LE.ZERO )
221: $ NEGCNT = NEGCNT + 1
222: *
223: DO 20 I = 2, N
224: TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
225: IF( ABS( TMP1 ).LT.PIVMIN )
226: $ TMP1 = -PIVMIN
227: IF( TMP1.LE.ZERO )
228: $ NEGCNT = NEGCNT + 1
229: 20 CONTINUE
230:
231: IF(NEGCNT.GE.IW) THEN
232: RIGHT = MID
233: ELSE
234: LEFT = MID
235: ENDIF
236: GOTO 10
237:
238: 30 CONTINUE
239: *
240: * Converged or maximum number of iterations reached
241: *
242: W = HALF * (LEFT + RIGHT)
243: WERR = HALF * ABS( RIGHT - LEFT )
244:
245: RETURN
246: *
247: * End of DLARRK
248: *
249: END
CVSweb interface <joel.bertrand@systella.fr>