LAPACK  3.6.1
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 September 2012
154 *
155 *> \ingroup auxOTHERauxiliary
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.4.2) --
172 * -- LAPACK is a software package provided by Univ. of Tennessee, --
173 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174 * September 2012
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  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
208  $ log( two ) ) + 2
209 *
210 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
211 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
212 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
213 * for an unconverged interval is set to the index of the next unconverged
214 * interval, and is -1 or 0 for a converged interval. Thus a linked
215 * list of unconverged intervals is set up.
216 *
217 
218  i1 = ifirst
219  i2 = ilast
220 * The number of unconverged intervals
221  nint = 0
222 * The last unconverged interval found
223  prev = 0
224  DO 75 i = i1, i2
225  k = 2*i
226  ii = i - offset
227  left = w( ii ) - werr( ii )
228  mid = w(ii)
229  right = w( ii ) + werr( ii )
230  width = right - mid
231  tmp = max( abs( left ), abs( right ) )
232 
233 * The following test prevents the test of converged intervals
234  IF( width.LT.rtol*tmp ) THEN
235 * This interval has already converged and does not need refinement.
236 * (Note that the gaps might change through refining the
237 * eigenvalues, however, they can only get bigger.)
238 * Remove it from the list.
239  iwork( k-1 ) = -1
240 * Make sure that I1 always points to the first unconverged interval
241  IF((i.EQ.i1).AND.(i.LT.i2)) i1 = i + 1
242  IF((prev.GE.i1).AND.(i.LE.i2)) iwork( 2*prev-1 ) = i + 1
243  ELSE
244 * unconverged interval found
245  prev = i
246 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
247 *
248 * Do while( CNT(LEFT).GT.I-1 )
249 *
250  fac = one
251  20 CONTINUE
252  cnt = 0
253  s = left
254  dplus = d( 1 ) - s
255  IF( dplus.LT.zero ) cnt = cnt + 1
256  DO 30 j = 2, n
257  dplus = d( j ) - s - e2( j-1 )/dplus
258  IF( dplus.LT.zero ) cnt = cnt + 1
259  30 CONTINUE
260  IF( cnt.GT.i-1 ) THEN
261  left = left - werr( ii )*fac
262  fac = two*fac
263  GO TO 20
264  END IF
265 *
266 * Do while( CNT(RIGHT).LT.I )
267 *
268  fac = one
269  50 CONTINUE
270  cnt = 0
271  s = right
272  dplus = d( 1 ) - s
273  IF( dplus.LT.zero ) cnt = cnt + 1
274  DO 60 j = 2, n
275  dplus = d( j ) - s - e2( j-1 )/dplus
276  IF( dplus.LT.zero ) cnt = cnt + 1
277  60 CONTINUE
278  IF( cnt.LT.i ) THEN
279  right = right + werr( ii )*fac
280  fac = two*fac
281  GO TO 50
282  END IF
283  nint = nint + 1
284  iwork( k-1 ) = i + 1
285  iwork( k ) = cnt
286  END IF
287  work( k-1 ) = left
288  work( k ) = right
289  75 CONTINUE
290 
291 
292  savi1 = i1
293 *
294 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
295 * and while (ITER.LT.MAXITR)
296 *
297  iter = 0
298  80 CONTINUE
299  prev = i1 - 1
300  i = i1
301  olnint = nint
302 
303  DO 100 p = 1, olnint
304  k = 2*i
305  ii = i - offset
306  next = iwork( k-1 )
307  left = work( k-1 )
308  right = work( k )
309  mid = half*( left + right )
310 
311 * semiwidth of interval
312  width = right - mid
313  tmp = max( abs( left ), abs( right ) )
314 
315  IF( ( width.LT.rtol*tmp ) .OR.
316  $ (iter.EQ.maxitr) )THEN
317 * reduce number of unconverged intervals
318  nint = nint - 1
319 * Mark interval as converged.
320  iwork( k-1 ) = 0
321  IF( i1.EQ.i ) THEN
322  i1 = next
323  ELSE
324 * Prev holds the last unconverged interval previously examined
325  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
326  END IF
327  i = next
328  GO TO 100
329  END IF
330  prev = i
331 *
332 * Perform one bisection step
333 *
334  cnt = 0
335  s = mid
336  dplus = d( 1 ) - s
337  IF( dplus.LT.zero ) cnt = cnt + 1
338  DO 90 j = 2, n
339  dplus = d( j ) - s - e2( j-1 )/dplus
340  IF( dplus.LT.zero ) cnt = cnt + 1
341  90 CONTINUE
342  IF( cnt.LE.i-1 ) THEN
343  work( k-1 ) = mid
344  ELSE
345  work( k ) = mid
346  END IF
347  i = next
348 
349  100 CONTINUE
350  iter = iter + 1
351 * do another loop if there are still unconverged intervals
352 * However, in the last iteration, all intervals are accepted
353 * since this is the best we can do.
354  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
355 *
356 *
357 * At this point, all the intervals have converged
358  DO 110 i = savi1, ilast
359  k = 2*i
360  ii = i - offset
361 * All intervals marked by '0' have been refined.
362  IF( iwork( k-1 ).EQ.0 ) THEN
363  w( ii ) = half*( work( k-1 )+work( k ) )
364  werr( ii ) = work( k ) - w( ii )
365  END IF
366  110 CONTINUE
367 *
368 
369  RETURN
370 *
371 * End of DLARRJ
372 *
373  END
dlarrj
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