1DLAEBZ(1) LAPACK auxiliary routine (version 3.2) DLAEBZ(1)
2
6 DLAEBZ - contains the iteration loops which compute and use the func‐
7 tion N(w), which is the count of eigenvalues of a symmetric tridiagonal
8 matrix T less than or equal to its argument w
9
11 SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL,
12 PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK,
13 IWORK, INFO )
14 INTEGER IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
16 DOUBLE PRECISION ABSTOL, PIVMIN, RELTOL
18 INTEGER IWORK( * ), NAB( MMAX, * ), NVAL( * )
20 DOUBLE PRECISION AB( MMAX, * ), C( * ), D( * ), E( * ), E2(
22 * ), WORK( * )
23
25 DLAEBZ contains the iteration loops which compute and use the function
26 N(w), which is the count of eigenvalues of a symmetric tridiagonal
27 matrix T less than or equal to its argument w. It performs a choice
28 of two types of loops:
29 IJOB=1, followed by
30 IJOB=2: It takes as input a list of intervals and returns a list of
31 sufficiently small intervals whose union contains the same
32 eigenvalues as the union of the original intervals.
33 The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
34 The output interval (AB(j,1),AB(j,2)] will contain
35 eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
36 IJOB=3: It performs a binary search in each input interval
37 (AB(j,1),AB(j,2)] for a point w(j) such that
38 N(w(j))=NVAL(j), and uses C(j) as the starting point of
39 the search. If such a w(j) is found, then on output
40 AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output
41 (AB(j,1),AB(j,2)] will be a small interval containing the
42 point where N(w) jumps through NVAL(j), unless that point
43 lies outside the initial interval.
44 Note that the intervals are in all cases half-open intervals, i.e., of
45 the form (a,b] , which includes b but not a . To avoid underflow,
46 the matrix should be scaled so that its largest element is no greater
47 than overflow**(1/2) * underflow**(1/4) in absolute value. To assure
48 the most accurate computation of small eigenvalues, the matrix should
49 be scaled to be
50 not much smaller than that, either.
51 See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal Matrix",
52 Report CS41, Computer Science Dept., Stanford
53 University, July 21, 1966
54 Note: the arguments are, in general, *not* checked for unreasonable
55 values.
56
58 IJOB (input) INTEGER
59 Specifies what is to be done:
60 = 1: Compute NAB for the initial intervals.
61 = 2: Perform bisection iteration to find eigenvalues of T.
62 = 3: Perform bisection iteration to invert N(w), i.e., to find
63 a point which has a specified number of eigenvalues of T to its
64 left. Other values will cause DLAEBZ to return with INFO=-1.
65 NITMAX (input) INTEGER
67 The maximum number of "levels" of bisection to be performed,
68 i.e., an interval of width W will not be made smaller than
69 2^(-NITMAX) * W. If not all intervals have converged after
70 NITMAX iterations, then INFO is set to the number of non-con‐
71 verged intervals.
72 N (input) INTEGER
74 The dimension n of the tridiagonal matrix T. It must be at
75 least 1.
76 MMAX (input) INTEGER
78 The maximum number of intervals. If more than MMAX intervals
79 are generated, then DLAEBZ will quit with INFO=MMAX+1.
80 MINP (input) INTEGER
82 The initial number of intervals. It may not be greater than
83 MMAX.
84 NBMIN (input) INTEGER
86 The smallest number of intervals that should be processed using
87 a vector loop. If zero, then only the scalar loop will be
88 used.
89 ABSTOL (input) DOUBLE PRECISION
91 The minimum (absolute) width of an interval. When an interval
92 is narrower than ABSTOL, or than RELTOL times the larger (in
93 magnitude) endpoint, then it is considered to be sufficiently
94 small, i.e., converged. This must be at least zero.
95 RELTOL (input) DOUBLE PRECISION
97 The minimum relative width of an interval. When an interval is
98 narrower than ABSTOL, or than RELTOL times the larger (in mag‐
99 nitude) endpoint, then it is considered to be sufficiently
100 small, i.e., converged. Note: this should always be at least
101 radix*machine epsilon.
102 PIVMIN (input) DOUBLE PRECISION
104 The minimum absolute value of a "pivot" in the Sturm sequence
105 loop. This *must* be at least max |e(j)**2| * safe_min and
106 at least safe_min, where safe_min is at least the smallest num‐
107 ber that can divide one without overflow.
108 D (input) DOUBLE PRECISION array, dimension (N)
110 The diagonal elements of the tridiagonal matrix T.
111 E (input) DOUBLE PRECISION array, dimension (N)
113 The offdiagonal elements of the tridiagonal matrix T in posi‐
114 tions 1 through N-1. E(N) is arbitrary.
115 E2 (input) DOUBLE PRECISION array, dimension (N)
117 The squares of the offdiagonal elements of the tridiagonal
118 matrix T. E2(N) is ignored.
119 NVAL (input/output) INTEGER array, dimension (MINP)
121 If IJOB=1 or 2, not referenced. If IJOB=3, the desired values
122 of N(w). The elements of NVAL will be reordered to correspond
123 with the intervals in AB. Thus, NVAL(j) on output will not, in
124 general be the same as NVAL(j) on input, but it will correspond
125 with the interval (AB(j,1),AB(j,2)] on output.
126 AB (input/output) DOUBLE PRECISION array, dimension (MMAX,2)
128 The endpoints of the intervals. AB(j,1) is a(j), the left
129 endpoint of the j-th interval, and AB(j,2) is b(j), the right
130 endpoint of the j-th interval. The input intervals will, in
131 general, be modified, split, and reordered by the calculation.
132 C (input/output) DOUBLE PRECISION array, dimension (MMAX)
134 If IJOB=1, ignored. If IJOB=2, workspace. If IJOB=3, then on
135 input C(j) should be initialized to the first search point in
136 the binary search.
137 MOUT (output) INTEGER
139 If IJOB=1, the number of eigenvalues in the intervals. If
140 IJOB=2 or 3, the number of intervals output. If IJOB=3, MOUT
141 will equal MINP.
142 NAB (input/output) INTEGER array, dimension (MMAX,2)
144 If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
145 If IJOB=2, then on input, NAB(i,j) should be set. It must sat‐
146 isfy the condition: N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <=
147 N(AB(i,2)), which means that in interval i only eigenvalues
148 NAB(i,1)+1,...,NAB(i,2) will be considered. Usually,
149 NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
150 IJOB=1. On output, NAB(i,j) will contain
151 max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of the
152 input interval that the output interval (AB(j,1),AB(j,2)] came
153 from, and na(k) and nb(k) are the the input values of NAB(k,1)
154 and NAB(k,2). If IJOB=3, then on output, NAB(i,j) contains
155 N(AB(i,j)), unless N(w) > NVAL(i) for all search points w , in
156 which case NAB(i,1) will not be modified, i.e., the output
157 value will be the same as the input value (modulo reorderings
158 -- see NVAL and AB), or unless N(w) < NVAL(i) for all search
159 points w , in which case NAB(i,2) will not be modified. Nor‐
160 mally, NAB should be set to some distinctive value(s) before
161 DLAEBZ is called.
162 WORK (workspace) DOUBLE PRECISION array, dimension (MMAX)
164 Workspace.
165 IWORK (workspace) INTEGER array, dimension (MMAX)
167 Workspace.
168 INFO (output) INTEGER
170 = 0: All intervals converged.
171 = 1--MMAX: The last INFO intervals did not converge.
172 = MMAX+1: More than MMAX intervals were generated.
173
175 This routine is intended to be called only by other LAPACK rou‐
176 tines, thus the interface is less user-friendly. It is intended for
177 two purposes:
178 (a) finding eigenvalues. In this case, DLAEBZ should have one or
179 more initial intervals set up in AB, and DLAEBZ should be called
180 with IJOB=1. This sets up NAB, and also counts the eigenvalues.
181 Intervals with no eigenvalues would usually be thrown out at
182 this point. Also, if not all the eigenvalues in an interval i
183 are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
184 For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
185 eigenvalue. DLAEBZ is then called with IJOB=2 and MMAX
186 no smaller than the value of MOUT returned by the call with
187 IJOB=1. After this (IJOB=2) call, eigenvalues NAB(i,1)+1
188 through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
189 tolerance specified by ABSTOL and RELTOL.
190 (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
191 In this case, start with a Gershgorin interval (a,b). Set up
192 AB to contain 2 search intervals, both initially (a,b). One
193 NVAL element should contain f-1 and the other should contain l
194 , while C should contain a and b, resp. NAB(i,1) should be -1
195 and NAB(i,2) should be N+1, to flag an error if the desired
196 interval does not lie in (a,b). DLAEBZ is then called with
197 IJOB=3. On exit, if w(f-1) < w(f), then one of the intervals --
198 j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
199 if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
200 >= 0, then the interval will have N(AB(j,1))=NAB(j,1)=f-k and
201 N(AB(j,2))=NAB(j,2)=f+r. The cases w(l) < w(l+1) and
202 w(l-r)=...=w(l+k) are handled similarly.
203 LAPACK auxiliary routine (versionNo3v.e2m)ber 2008 DLAEBZ(1)