LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlarrb.f
Go to the documentation of this file.
1 *> \brief \b DLARRB 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 DLARRB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrb.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrb.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrb.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRB( 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 * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * DOUBLE PRECISION 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, DLARRB
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
87 *> \endverbatim
88 *>
89 *> \param[in] RTOL2
90 *> \verbatim
91 *> RTOL2 is DOUBLE PRECISION
92 *> Tolerance for the convergence of the bisection intervals.
93 *> An interval [LEFT,RIGHT] has converged if
94 *> RIGHT-LEFT.LT.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 DOUBLE PRECISION 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 throug
111 *> ILAST.
112 *> On output, these estimates are refined.
113 *> \endverbatim
114 *>
115 *> \param[in,out] WGAP
116 *> \verbatim
117 *> WGAP is DOUBLE PRECISION 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.EQ.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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
148 *> The minimum pivot in the Sturm sequence.
149 *> \endverbatim
150 *>
151 *> \param[in] SPDIAM
152 *> \verbatim
153 *> SPDIAM is DOUBLE PRECISION
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 September 2012
182 *
183 *> \ingroup auxOTHERauxiliary
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 dlarrb( 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.4.2) --
200 * -- LAPACK is a software package provided by Univ. of Tennessee, --
201 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
202 * September 2012
203 *
204 * .. Scalar Arguments ..
205  INTEGER ifirst, ilast, info, n, offset, twist
206  DOUBLE PRECISION pivmin, rtol1, rtol2, spdiam
207 * ..
208 * .. Array Arguments ..
209  INTEGER iwork( * )
210  DOUBLE PRECISION d( * ), lld( * ), w( * ),
211  $ werr( * ), wgap( * ), work( * )
212 * ..
213 *
214 * =====================================================================
215 *
216 * .. Parameters ..
217  DOUBLE PRECISION zero, two, half
218  parameter( zero = 0.0d0, two = 2.0d0,
219  $ half = 0.5d0 )
220  INTEGER maxitr
221 * ..
222 * .. Local Scalars ..
223  INTEGER i, i1, ii, ip, iter, k, negcnt, next, nint,
224  $ olnint, prev, r
225  DOUBLE PRECISION back, cvrgd, gap, left, lgap, mid, mnwdth,
226  $ rgap, right, tmp, width
227 * ..
228 * .. External Functions ..
229  INTEGER dlaneg
230  EXTERNAL dlaneg
231 *
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC abs, max, min
235 * ..
236 * .. Executable Statements ..
237 *
238  info = 0
239 *
240  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
241  $ log( two ) ) + 2
242  mnwdth = two * pivmin
243 *
244  r = twist
245  IF((r.LT.1).OR.(r.GT.n)) r = n
246 *
247 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
248 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
249 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
250 * for an unconverged interval is set to the index of the next unconverged
251 * interval, and is -1 or 0 for a converged interval. Thus a linked
252 * list of unconverged intervals is set up.
253 *
254  i1 = ifirst
255 * The number of unconverged intervals
256  nint = 0
257 * The last unconverged interval found
258  prev = 0
259 
260  rgap = wgap( i1-offset )
261  DO 75 i = i1, ilast
262  k = 2*i
263  ii = i - offset
264  left = w( ii ) - werr( ii )
265  right = w( ii ) + werr( ii )
266  lgap = rgap
267  rgap = wgap( ii )
268  gap = min( lgap, rgap )
269 
270 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
271 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
272 *
273 * Do while( NEGCNT(LEFT).GT.I-1 )
274 *
275  back = werr( ii )
276  20 continue
277  negcnt = dlaneg( n, d, lld, left, pivmin, r )
278  IF( negcnt.GT.i-1 ) THEN
279  left = left - back
280  back = two*back
281  go to 20
282  END IF
283 *
284 * Do while( NEGCNT(RIGHT).LT.I )
285 * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
286 *
287  back = werr( ii )
288  50 continue
289 
290  negcnt = dlaneg( n, d, lld, right, pivmin, r )
291  IF( negcnt.LT.i ) THEN
292  right = right + back
293  back = two*back
294  go to 50
295  END IF
296  width = half*abs( left - right )
297  tmp = max( abs( left ), abs( right ) )
298  cvrgd = max(rtol1*gap,rtol2*tmp)
299  IF( width.LE.cvrgd .OR. width.LE.mnwdth ) THEN
300 * This interval has already converged and does not need refinement.
301 * (Note that the gaps might change through refining the
302 * eigenvalues, however, they can only get bigger.)
303 * Remove it from the list.
304  iwork( k-1 ) = -1
305 * Make sure that I1 always points to the first unconverged interval
306  IF((i.EQ.i1).AND.(i.LT.ilast)) i1 = i + 1
307  IF((prev.GE.i1).AND.(i.LE.ilast)) iwork( 2*prev-1 ) = i + 1
308  ELSE
309 * unconverged interval found
310  prev = i
311  nint = nint + 1
312  iwork( k-1 ) = i + 1
313  iwork( k ) = negcnt
314  END IF
315  work( k-1 ) = left
316  work( k ) = right
317  75 continue
318 
319 *
320 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
321 * and while (ITER.LT.MAXITR)
322 *
323  iter = 0
324  80 continue
325  prev = i1 - 1
326  i = i1
327  olnint = nint
328 
329  DO 100 ip = 1, olnint
330  k = 2*i
331  ii = i - offset
332  rgap = wgap( ii )
333  lgap = rgap
334  IF(ii.GT.1) lgap = wgap( ii-1 )
335  gap = min( lgap, rgap )
336  next = iwork( k-1 )
337  left = work( k-1 )
338  right = work( k )
339  mid = half*( left + right )
340 
341 * semiwidth of interval
342  width = right - mid
343  tmp = max( abs( left ), abs( right ) )
344  cvrgd = max(rtol1*gap,rtol2*tmp)
345  IF( ( width.LE.cvrgd ) .OR. ( width.LE.mnwdth ).OR.
346  $ ( iter.EQ.maxitr ) )THEN
347 * reduce number of unconverged intervals
348  nint = nint - 1
349 * Mark interval as converged.
350  iwork( k-1 ) = 0
351  IF( i1.EQ.i ) THEN
352  i1 = next
353  ELSE
354 * Prev holds the last unconverged interval previously examined
355  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
356  END IF
357  i = next
358  go to 100
359  END IF
360  prev = i
361 *
362 * Perform one bisection step
363 *
364  negcnt = dlaneg( n, d, lld, mid, pivmin, r )
365  IF( negcnt.LE.i-1 ) THEN
366  work( k-1 ) = mid
367  ELSE
368  work( k ) = mid
369  END IF
370  i = next
371  100 continue
372  iter = iter + 1
373 * do another loop if there are still unconverged intervals
374 * However, in the last iteration, all intervals are accepted
375 * since this is the best we can do.
376  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) go to 80
377 *
378 *
379 * At this point, all the intervals have converged
380  DO 110 i = ifirst, ilast
381  k = 2*i
382  ii = i - offset
383 * All intervals marked by '0' have been refined.
384  IF( iwork( k-1 ).EQ.0 ) THEN
385  w( ii ) = half*( work( k-1 )+work( k ) )
386  werr( ii ) = work( k ) - w( ii )
387  END IF
388  110 continue
389 *
390  DO 111 i = ifirst+1, ilast
391  k = 2*i
392  ii = i - offset
393  wgap( ii-1 ) = max( zero,
394  $ w(ii) - werr(ii) - w( ii-1 ) - werr( ii-1 ))
395  111 continue
396 
397  return
398 *
399 * End of DLARRB
400 *
401  END