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