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