LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slarrj.f
Go to the documentation of this file.
1 *> \brief \b SLARRJ 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 SLARRJ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrj.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrj.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrj.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARRJ( 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 * REAL PIVMIN, RTOL, SPDIAM
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * REAL 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, SLARRJ
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 REAL array, dimension (N)
62 *> The N diagonal elements of T.
63 *> \endverbatim
64 *>
65 *> \param[in] E2
66 *> \verbatim
67 *> E2 is REAL 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 REAL
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 REAL 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 REAL 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 REAL 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 REAL
130 *> The minimum pivot in the Sturm sequence for T.
131 *> \endverbatim
132 *>
133 *> \param[in] SPDIAM
134 *> \verbatim
135 *> SPDIAM is REAL
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 slarrj( 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  REAL pivmin, rtol, spdiam
179 * ..
180 * .. Array Arguments ..
181  INTEGER iwork( * )
182  REAL d( * ), e2( * ), w( * ),
183  $ werr( * ), work( * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  REAL zero, one, two, half
190  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
191  $ half = 0.5e0 )
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  REAL 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 SLARRJ
372 *
373  END