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: * =====================================================================
145: SUBROUTINE DLARRK( N, IW, GL, GU,
146: $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
147: *
148: * -- LAPACK auxiliary routine (version 3.4.0) --
149: * -- LAPACK is a software package provided by Univ. of Tennessee, --
150: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
151: * November 2011
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>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrk.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack