LAPACK  3.9.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 *> \date June 2017
182 *
183 *> \ingroup OTHERauxiliary
184 *
185 *> \par Contributors:
186 * ==================
187 *>
188 *> Beresford Parlett, University of California, Berkeley, USA \n
189 *> Jim Demmel, University of California, Berkeley, USA \n
190 *> Inderjit Dhillon, University of Texas, Austin, USA \n
191 *> Osni Marques, LBNL/NERSC, USA \n
192 *> Christof Voemel, University of California, Berkeley, USA
193 *
194 * =====================================================================
195  SUBROUTINE slarrb( N, D, LLD, IFIRST, ILAST, RTOL1,
196  $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
197  $ PIVMIN, SPDIAM, TWIST, INFO )
198 *
199 * -- LAPACK auxiliary routine (version 3.7.1) --
200 * -- LAPACK is a software package provided by Univ. of Tennessee, --
201 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
202 * June 2017
203 *
204 * .. Scalar Arguments ..
205  INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
206  REAL PIVMIN, RTOL1, RTOL2, SPDIAM
207 * ..
208 * .. Array Arguments ..
209  INTEGER IWORK( * )
210  REAL D( * ), LLD( * ), W( * ),
211  $ werr( * ), wgap( * ), work( * )
212 * ..
213 *
214 * =====================================================================
215 *
216 * .. Parameters ..
217  REAL ZERO, TWO, HALF
218  PARAMETER ( ZERO = 0.0e0, two = 2.0e0,
219  $ half = 0.5e0 )
220  INTEGER MAXITR
221 * ..
222 * .. Local Scalars ..
223  INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
224  $ OLNINT, PREV, R
225  REAL BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
226  $ RGAP, RIGHT, TMP, WIDTH
227 * ..
228 * .. External Functions ..
229  INTEGER SLANEG
230  EXTERNAL SLANEG
231 *
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, max, min
235 * ..
236 * .. Executable Statements ..
237 *
238  info = 0
239 *
240 * Quick return if possible
241 *
242  IF( n.LE.0 ) THEN
243  RETURN
244  END IF
245 *
246  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
247  $ log( two ) ) + 2
248  mnwdth = two * pivmin
249 *
250  r = twist
251  IF((r.LT.1).OR.(r.GT.n)) r = n
252 *
253 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
254 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
255 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
256 * for an unconverged interval is set to the index of the next unconverged
257 * interval, and is -1 or 0 for a converged interval. Thus a linked
258 * list of unconverged intervals is set up.
259 *
260  i1 = ifirst
261 * The number of unconverged intervals
262  nint = 0
263 * The last unconverged interval found
264  prev = 0
265 
266  rgap = wgap( i1-offset )
267  DO 75 i = i1, ilast
268  k = 2*i
269  ii = i - offset
270  left = w( ii ) - werr( ii )
271  right = w( ii ) + werr( ii )
272  lgap = rgap
273  rgap = wgap( ii )
274  gap = min( lgap, rgap )
275 
276 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
277 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
278 *
279 * Do while( NEGCNT(LEFT).GT.I-1 )
280 *
281  back = werr( ii )
282  20 CONTINUE
283  negcnt = slaneg( n, d, lld, left, pivmin, r )
284  IF( negcnt.GT.i-1 ) THEN
285  left = left - back
286  back = two*back
287  GO TO 20
288  END IF
289 *
290 * Do while( NEGCNT(RIGHT).LT.I )
291 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
292 *
293  back = werr( ii )
294  50 CONTINUE
295 
296  negcnt = slaneg( n, d, lld, right, pivmin, r )
297  IF( negcnt.LT.i ) THEN
298  right = right + back
299  back = two*back
300  GO TO 50
301  END IF
302  width = half*abs( left - right )
303  tmp = max( abs( left ), abs( right ) )
304  cvrgd = max(rtol1*gap,rtol2*tmp)
305  IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
306 * This interval has already converged and does not need refinement.
307 * (Note that the gaps might change through refining the
308 * eigenvalues, however, they can only get bigger.)
309 * Remove it from the list.
310  iwork( k-1 ) = -1
311 * Make sure that I1 always points to the first unconverged interval
312  IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
313  IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
314  ELSE
315 * unconverged interval found
316  prev = i
317  nint = nint + 1
318  iwork( k-1 ) = i + 1
319  iwork( k ) = negcnt
320  END IF
321  work( k-1 ) = left
322  work( k ) = right
323  75 CONTINUE
324 
325 *
326 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
327 * and while (ITER.LT.MAXITR)
328 *
329  iter = 0
330  80 CONTINUE
331  prev = i1 - 1
332  i = i1
333  olnint = nint
334 
335  DO 100 ip = 1, olnint
336  k = 2*i
337  ii = i - offset
338  rgap = wgap( ii )
339  lgap = rgap
340  IF(ii.GT.1) lgap = wgap( ii-1 )
341  gap = min( lgap, rgap )
342  next = iwork( k-1 )
343  left = work( k-1 )
344  right = work( k )
345  mid = half*( left + right )
346 
347 * semiwidth of interval
348  width = right - mid
349  tmp = max( abs( left ), abs( right ) )
350  cvrgd = max(rtol1*gap,rtol2*tmp)
351  IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
352  $ ( iter.EQ.maxitr ) )THEN
353 * reduce number of unconverged intervals
354  nint = nint - 1
355 * Mark interval as converged.
356  iwork( k-1 ) = 0
357  IF( i1.EQ.i ) THEN
358  i1 = next
359  ELSE
360 * Prev holds the last unconverged interval previously examined
361  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
362  END IF
363  i = next
364  GO TO 100
365  END IF
366  prev = i
367 *
368 * Perform one bisection step
369 *
370  negcnt = slaneg( n, d, lld, mid, pivmin, r )
371  IF( negcnt.LE.i-1 ) THEN
372  work( k-1 ) = mid
373  ELSE
374  work( k ) = mid
375  END IF
376  i = next
377  100 CONTINUE
378  iter = iter + 1
379 * do another loop if there are still unconverged intervals
380 * However, in the last iteration, all intervals are accepted
381 * since this is the best we can do.
382  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
383 *
384 *
385 * At this point, all the intervals have converged
386  DO 110 i = ifirst, ilast
387  k = 2*i
388  ii = i - offset
389 * All intervals marked by '0' have been refined.
390  IF( iwork( k-1 ).EQ.0 ) THEN
391  w( ii ) = half*( work( k-1 )+work( k ) )
392  werr( ii ) = work( k ) - w( ii )
393  END IF
394  110 CONTINUE
395 *
396  DO 111 i = ifirst+1, ilast
397  k = 2*i
398  ii = i - offset
399  wgap( ii-1 ) = max( zero,
400  $ w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
401  111 CONTINUE
402 
403  RETURN
404 *
405 * End of SLARRB
406 *
407  END
slarrb
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:198