LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slarrb.f
Go to the documentation of this file.
1 *> \brief \b SLARRB provides limited bisection to locate eigenvalues for more accuracy.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLARRB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrb.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrb.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrb.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
22 * RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
23 * PIVMIN, SPDIAM, TWIST, INFO )
24 *
25 * .. Scalar Arguments ..
26 * INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
27 * REAL PIVMIN, RTOL1, RTOL2, SPDIAM
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * REAL D( * ), LLD( * ), W( * ),
32 * $ WERR( * ), WGAP( * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> Given the relatively robust representation(RRR) L D L^T, SLARRB
42 *> does "limited" bisection to refine the eigenvalues of L D L^T,
43 *> W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
44 *> guesses for these eigenvalues are input in W, the corresponding estimate
45 *> of the error in these guesses and their gaps are input in WERR
46 *> and WGAP, respectively. During bisection, intervals
47 *> [left, right] are maintained by storing their mid-points and
48 *> semi-widths in the arrays W and WERR respectively.
49 *> \endverbatim
50 *
51 * Arguments:
52 * ==========
53 *
54 *> \param[in] N
55 *> \verbatim
56 *> N is INTEGER
57 *> The order of the matrix.
58 *> \endverbatim
59 *>
60 *> \param[in] D
61 *> \verbatim
62 *> D is REAL array, dimension (N)
63 *> The N diagonal elements of the diagonal matrix D.
64 *> \endverbatim
65 *>
66 *> \param[in] LLD
67 *> \verbatim
68 *> LLD is REAL array, dimension (N-1)
69 *> The (N-1) elements L(i)*L(i)*D(i).
70 *> \endverbatim
71 *>
72 *> \param[in] IFIRST
73 *> \verbatim
74 *> IFIRST is INTEGER
75 *> The index of the first eigenvalue to be computed.
76 *> \endverbatim
77 *>
78 *> \param[in] ILAST
79 *> \verbatim
80 *> ILAST is INTEGER
81 *> The index of the last eigenvalue to be computed.
82 *> \endverbatim
83 *>
84 *> \param[in] RTOL1
85 *> \verbatim
86 *> RTOL1 is REAL
87 *> \endverbatim
88 *>
89 *> \param[in] RTOL2
90 *> \verbatim
91 *> RTOL2 is REAL
92 *> Tolerance for the convergence of the bisection intervals.
93 *> An interval [LEFT,RIGHT] has converged if
94 *> RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
95 *> where GAP is the (estimated) distance to the nearest
96 *> eigenvalue.
97 *> \endverbatim
98 *>
99 *> \param[in] OFFSET
100 *> \verbatim
101 *> OFFSET is INTEGER
102 *> Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
103 *> through ILAST-OFFSET elements of these arrays are to be used.
104 *> \endverbatim
105 *>
106 *> \param[in,out] W
107 *> \verbatim
108 *> W is REAL array, dimension (N)
109 *> On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
110 *> estimates of the eigenvalues of L D L^T indexed IFIRST through
111 *> ILAST.
112 *> On output, these estimates are refined.
113 *> \endverbatim
114 *>
115 *> \param[in,out] WGAP
116 *> \verbatim
117 *> WGAP is REAL array, dimension (N-1)
118 *> On input, the (estimated) gaps between consecutive
119 *> eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
120 *> eigenvalues I and I+1. Note that if IFIRST = ILAST
121 *> then WGAP(IFIRST-OFFSET) must be set to ZERO.
122 *> On output, these gaps are refined.
123 *> \endverbatim
124 *>
125 *> \param[in,out] WERR
126 *> \verbatim
127 *> WERR is REAL array, dimension (N)
128 *> On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
129 *> the errors in the estimates of the corresponding elements in W.
130 *> On output, these errors are refined.
131 *> \endverbatim
132 *>
133 *> \param[out] WORK
134 *> \verbatim
135 *> WORK is REAL array, dimension (2*N)
136 *> Workspace.
137 *> \endverbatim
138 *>
139 *> \param[out] IWORK
140 *> \verbatim
141 *> IWORK is INTEGER array, dimension (2*N)
142 *> Workspace.
143 *> \endverbatim
144 *>
145 *> \param[in] PIVMIN
146 *> \verbatim
147 *> PIVMIN is REAL
148 *> The minimum pivot in the Sturm sequence.
149 *> \endverbatim
150 *>
151 *> \param[in] SPDIAM
152 *> \verbatim
153 *> SPDIAM is REAL
154 *> The spectral diameter of the matrix.
155 *> \endverbatim
156 *>
157 *> \param[in] TWIST
158 *> \verbatim
159 *> TWIST is INTEGER
160 *> The twist index for the twisted factorization that is used
161 *> for the negcount.
162 *> TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
163 *> TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
164 *> TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
165 *> \endverbatim
166 *>
167 *> \param[out] INFO
168 *> \verbatim
169 *> INFO is INTEGER
170 *> Error flag.
171 *> \endverbatim
172 *
173 * Authors:
174 * ========
175 *
176 *> \author Univ. of Tennessee
177 *> \author Univ. of California Berkeley
178 *> \author Univ. of Colorado Denver
179 *> \author NAG Ltd.
180 *
181 *> \ingroup OTHERauxiliary
182 *
183 *> \par Contributors:
184 * ==================
185 *>
186 *> Beresford Parlett, University of California, Berkeley, USA \n
187 *> Jim Demmel, University of California, Berkeley, USA \n
188 *> Inderjit Dhillon, University of Texas, Austin, USA \n
189 *> Osni Marques, LBNL/NERSC, USA \n
190 *> Christof Voemel, University of California, Berkeley, USA
191 *
192 * =====================================================================
193  SUBROUTINE slarrb( N, D, LLD, IFIRST, ILAST, RTOL1,
194  $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
195  $ PIVMIN, SPDIAM, TWIST, INFO )
196 *
197 * -- LAPACK auxiliary routine --
198 * -- LAPACK is a software package provided by Univ. of Tennessee, --
199 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
200 *
201 * .. Scalar Arguments ..
202  INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
203  REAL PIVMIN, RTOL1, RTOL2, SPDIAM
204 * ..
205 * .. Array Arguments ..
206  INTEGER IWORK( * )
207  REAL D( * ), LLD( * ), W( * ),
208  $ werr( * ), wgap( * ), work( * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  REAL ZERO, TWO, HALF
215  PARAMETER ( ZERO = 0.0e0, two = 2.0e0,
216  $ half = 0.5e0 )
217  INTEGER MAXITR
218 * ..
219 * .. Local Scalars ..
220  INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
221  $ OLNINT, PREV, R
222  REAL BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
223  $ RGAP, RIGHT, TMP, WIDTH
224 * ..
225 * .. External Functions ..
226  INTEGER SLANEG
227  EXTERNAL SLANEG
228 *
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC abs, max, min
232 * ..
233 * .. Executable Statements ..
234 *
235  info = 0
236 *
237 * Quick return if possible
238 *
239  IF( n.LE.0 ) THEN
240  RETURN
241  END IF
242 *
243  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
244  $ log( two ) ) + 2
245  mnwdth = two * pivmin
246 *
247  r = twist
248  IF((r.LT.1).OR.(r.GT.n)) r = n
249 *
250 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
251 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
252 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
253 * for an unconverged interval is set to the index of the next unconverged
254 * interval, and is -1 or 0 for a converged interval. Thus a linked
255 * list of unconverged intervals is set up.
256 *
257  i1 = ifirst
258 * The number of unconverged intervals
259  nint = 0
260 * The last unconverged interval found
261  prev = 0
262 
263  rgap = wgap( i1-offset )
264  DO 75 i = i1, ilast
265  k = 2*i
266  ii = i - offset
267  left = w( ii ) - werr( ii )
268  right = w( ii ) + werr( ii )
269  lgap = rgap
270  rgap = wgap( ii )
271  gap = min( lgap, rgap )
272 
273 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
274 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
275 *
276 * Do while( NEGCNT(LEFT).GT.I-1 )
277 *
278  back = werr( ii )
279  20 CONTINUE
280  negcnt = slaneg( n, d, lld, left, pivmin, r )
281  IF( negcnt.GT.i-1 ) THEN
282  left = left - back
283  back = two*back
284  GO TO 20
285  END IF
286 *
287 * Do while( NEGCNT(RIGHT).LT.I )
288 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
289 *
290  back = werr( ii )
291  50 CONTINUE
292 
293  negcnt = slaneg( n, d, lld, right, pivmin, r )
294  IF( negcnt.LT.i ) THEN
295  right = right + back
296  back = two*back
297  GO TO 50
298  END IF
299  width = half*abs( left - right )
300  tmp = max( abs( left ), abs( right ) )
301  cvrgd = max(rtol1*gap,rtol2*tmp)
302  IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
303 * This interval has already converged and does not need refinement.
304 * (Note that the gaps might change through refining the
305 * eigenvalues, however, they can only get bigger.)
306 * Remove it from the list.
307  iwork( k-1 ) = -1
308 * Make sure that I1 always points to the first unconverged interval
309  IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
310  IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
311  ELSE
312 * unconverged interval found
313  prev = i
314  nint = nint + 1
315  iwork( k-1 ) = i + 1
316  iwork( k ) = negcnt
317  END IF
318  work( k-1 ) = left
319  work( k ) = right
320  75 CONTINUE
321 
322 *
323 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
324 * and while (ITER.LT.MAXITR)
325 *
326  iter = 0
327  80 CONTINUE
328  prev = i1 - 1
329  i = i1
330  olnint = nint
331 
332  DO 100 ip = 1, olnint
333  k = 2*i
334  ii = i - offset
335  rgap = wgap( ii )
336  lgap = rgap
337  IF(ii.GT.1) lgap = wgap( ii-1 )
338  gap = min( lgap, rgap )
339  next = iwork( k-1 )
340  left = work( k-1 )
341  right = work( k )
342  mid = half*( left + right )
343 
344 * semiwidth of interval
345  width = right - mid
346  tmp = max( abs( left ), abs( right ) )
347  cvrgd = max(rtol1*gap,rtol2*tmp)
348  IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
349  $ ( iter.EQ.maxitr ) )THEN
350 * reduce number of unconverged intervals
351  nint = nint - 1
352 * Mark interval as converged.
353  iwork( k-1 ) = 0
354  IF( i1.EQ.i ) THEN
355  i1 = next
356  ELSE
357 * Prev holds the last unconverged interval previously examined
358  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
359  END IF
360  i = next
361  GO TO 100
362  END IF
363  prev = i
364 *
365 * Perform one bisection step
366 *
367  negcnt = slaneg( n, d, lld, mid, pivmin, r )
368  IF( negcnt.LE.i-1 ) THEN
369  work( k-1 ) = mid
370  ELSE
371  work( k ) = mid
372  END IF
373  i = next
374  100 CONTINUE
375  iter = iter + 1
376 * do another loop if there are still unconverged intervals
377 * However, in the last iteration, all intervals are accepted
378 * since this is the best we can do.
379  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
380 *
381 *
382 * At this point, all the intervals have converged
383  DO 110 i = ifirst, ilast
384  k = 2*i
385  ii = i - offset
386 * All intervals marked by '0' have been refined.
387  IF( iwork( k-1 ).EQ.0 ) THEN
388  w( ii ) = half*( work( k-1 )+work( k ) )
389  werr( ii ) = work( k ) - w( ii )
390  END IF
391  110 CONTINUE
392 *
393  DO 111 i = ifirst+1, ilast
394  k = 2*i
395  ii = i - offset
396  wgap( ii-1 ) = max( zero,
397  $ w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
398  111 CONTINUE
399 
400  RETURN
401 *
402 * End of SLARRB
403 *
404  END
subroutine slarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
SLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: slarrb.f:196