1: *> \brief \b DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
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: *> \ingroup OTHERauxiliary
141: *
142: * =====================================================================
143: SUBROUTINE DLARRK( N, IW, GL, GU,
144: $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
145: *
146: * -- LAPACK auxiliary routine --
147: * -- LAPACK is a software package provided by Univ. of Tennessee, --
148: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149: *
150: * .. Scalar Arguments ..
151: INTEGER INFO, IW, N
152: DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
153: * ..
154: * .. Array Arguments ..
155: DOUBLE PRECISION D( * ), E2( * )
156: * ..
157: *
158: * =====================================================================
159: *
160: * .. Parameters ..
161: DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
162: PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
163: $ FUDGE = TWO, ZERO = 0.0D0 )
164: * ..
165: * .. Local Scalars ..
166: INTEGER I, IT, ITMAX, NEGCNT
167: DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
168: $ TMP2, TNORM
169: * ..
170: * .. External Functions ..
171: DOUBLE PRECISION DLAMCH
172: EXTERNAL DLAMCH
173: * ..
174: * .. Intrinsic Functions ..
175: INTRINSIC ABS, INT, LOG, MAX
176: * ..
177: * .. Executable Statements ..
178: *
179: * Quick return if possible
180: *
181: IF( N.LE.0 ) THEN
182: INFO = 0
183: RETURN
184: END IF
185: *
186: * Get machine constants
187: EPS = DLAMCH( 'P' )
188:
189: TNORM = MAX( ABS( GL ), ABS( GU ) )
190: RTOLI = RELTOL
191: ATOLI = FUDGE*TWO*PIVMIN
192:
193: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
194: $ LOG( TWO ) ) + 2
195:
196: INFO = -1
197:
198: LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
199: RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
200: IT = 0
201:
202: 10 CONTINUE
203: *
204: * Check if interval converged or maximum number of iterations reached
205: *
206: TMP1 = ABS( RIGHT - LEFT )
207: TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
208: IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
209: INFO = 0
210: GOTO 30
211: ENDIF
212: IF(IT.GT.ITMAX)
213: $ GOTO 30
214:
215: *
216: * Count number of negative pivots for mid-point
217: *
218: IT = IT + 1
219: MID = HALF * (LEFT + RIGHT)
220: NEGCNT = 0
221: TMP1 = D( 1 ) - MID
222: IF( ABS( TMP1 ).LT.PIVMIN )
223: $ TMP1 = -PIVMIN
224: IF( TMP1.LE.ZERO )
225: $ NEGCNT = NEGCNT + 1
226: *
227: DO 20 I = 2, N
228: TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
229: IF( ABS( TMP1 ).LT.PIVMIN )
230: $ TMP1 = -PIVMIN
231: IF( TMP1.LE.ZERO )
232: $ NEGCNT = NEGCNT + 1
233: 20 CONTINUE
234:
235: IF(NEGCNT.GE.IW) THEN
236: RIGHT = MID
237: ELSE
238: LEFT = MID
239: ENDIF
240: GOTO 10
241:
242: 30 CONTINUE
243: *
244: * Converged or maximum number of iterations reached
245: *
246: W = HALF * (LEFT + RIGHT)
247: WERR = HALF * ABS( RIGHT - LEFT )
248:
249: RETURN
250: *
251: * End of DLARRK
252: *
253: 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