Annotation of rpl/lapack/lapack/dlarrk.f, revision 1.1.1.1
1.1 bertrand 1: SUBROUTINE DLARRK( N, IW, GL, GU,
2: $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
3: IMPLICIT NONE
4: *
5: * -- LAPACK auxiliary routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: INTEGER INFO, IW, N
12: DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION D( * ), E2( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLARRK computes one eigenvalue of a symmetric tridiagonal
22: * matrix T to suitable accuracy. This is an auxiliary code to be
23: * called from DSTEMR.
24: *
25: * To avoid overflow, the matrix must be scaled so that its
26: * largest element is no greater than overflow**(1/2) *
27: * underflow**(1/4) in absolute value, and for greatest
28: * accuracy, it should not be much smaller than that.
29: *
30: * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
31: * Matrix", Report CS41, Computer Science Dept., Stanford
32: * University, July 21, 1966.
33: *
34: * Arguments
35: * =========
36: *
37: * N (input) INTEGER
38: * The order of the tridiagonal matrix T. N >= 0.
39: *
40: * IW (input) INTEGER
41: * The index of the eigenvalues to be returned.
42: *
43: * GL (input) DOUBLE PRECISION
44: * GU (input) DOUBLE PRECISION
45: * An upper and a lower bound on the eigenvalue.
46: *
47: * D (input) DOUBLE PRECISION array, dimension (N)
48: * The n diagonal elements of the tridiagonal matrix T.
49: *
50: * E2 (input) DOUBLE PRECISION array, dimension (N-1)
51: * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
52: *
53: * PIVMIN (input) DOUBLE PRECISION
54: * The minimum pivot allowed in the Sturm sequence for T.
55: *
56: * RELTOL (input) DOUBLE PRECISION
57: * The minimum relative width of an interval. When an interval
58: * is narrower than RELTOL times the larger (in
59: * magnitude) endpoint, then it is considered to be
60: * sufficiently small, i.e., converged. Note: this should
61: * always be at least radix*machine epsilon.
62: *
63: * W (output) DOUBLE PRECISION
64: *
65: * WERR (output) DOUBLE PRECISION
66: * The error bound on the corresponding eigenvalue approximation
67: * in W.
68: *
69: * INFO (output) INTEGER
70: * = 0: Eigenvalue converged
71: * = -1: Eigenvalue did NOT converge
72: *
73: * Internal Parameters
74: * ===================
75: *
76: * FUDGE DOUBLE PRECISION, default = 2
77: * A "fudge factor" to widen the Gershgorin intervals.
78: *
79: * =====================================================================
80: *
81: * .. Parameters ..
82: DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
83: PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
84: $ FUDGE = TWO, ZERO = 0.0D0 )
85: * ..
86: * .. Local Scalars ..
87: INTEGER I, IT, ITMAX, NEGCNT
88: DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
89: $ TMP2, TNORM
90: * ..
91: * .. External Functions ..
92: DOUBLE PRECISION DLAMCH
93: EXTERNAL DLAMCH
94: * ..
95: * .. Intrinsic Functions ..
96: INTRINSIC ABS, INT, LOG, MAX
97: * ..
98: * .. Executable Statements ..
99: *
100: * Get machine constants
101: EPS = DLAMCH( 'P' )
102:
103: TNORM = MAX( ABS( GL ), ABS( GU ) )
104: RTOLI = RELTOL
105: ATOLI = FUDGE*TWO*PIVMIN
106:
107: ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
108: $ LOG( TWO ) ) + 2
109:
110: INFO = -1
111:
112: LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
113: RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
114: IT = 0
115:
116: 10 CONTINUE
117: *
118: * Check if interval converged or maximum number of iterations reached
119: *
120: TMP1 = ABS( RIGHT - LEFT )
121: TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
122: IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
123: INFO = 0
124: GOTO 30
125: ENDIF
126: IF(IT.GT.ITMAX)
127: $ GOTO 30
128:
129: *
130: * Count number of negative pivots for mid-point
131: *
132: IT = IT + 1
133: MID = HALF * (LEFT + RIGHT)
134: NEGCNT = 0
135: TMP1 = D( 1 ) - MID
136: IF( ABS( TMP1 ).LT.PIVMIN )
137: $ TMP1 = -PIVMIN
138: IF( TMP1.LE.ZERO )
139: $ NEGCNT = NEGCNT + 1
140: *
141: DO 20 I = 2, N
142: TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
143: IF( ABS( TMP1 ).LT.PIVMIN )
144: $ TMP1 = -PIVMIN
145: IF( TMP1.LE.ZERO )
146: $ NEGCNT = NEGCNT + 1
147: 20 CONTINUE
148:
149: IF(NEGCNT.GE.IW) THEN
150: RIGHT = MID
151: ELSE
152: LEFT = MID
153: ENDIF
154: GOTO 10
155:
156: 30 CONTINUE
157: *
158: * Converged or maximum number of iterations reached
159: *
160: W = HALF * (LEFT + RIGHT)
161: WERR = HALF * ABS( RIGHT - LEFT )
162:
163: RETURN
164: *
165: * End of DLARRK
166: *
167: 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