LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
dlarrb
subroutine dlarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: dlarrb.f:198