LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 < 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 *> \ingroup OTHERauxiliary
154 *
155 *> \par Contributors:
156 * ==================
157 *>
158 *> Beresford Parlett, University of California, Berkeley, USA \n
159 *> Jim Demmel, University of California, Berkeley, USA \n
160 *> Inderjit Dhillon, University of Texas, Austin, USA \n
161 *> Osni Marques, LBNL/NERSC, USA \n
162 *> Christof Voemel, University of California, Berkeley, USA
163 *
164 * =====================================================================
165  SUBROUTINE slarrj( N, D, E2, IFIRST, ILAST,
166  $ RTOL, OFFSET, W, WERR, WORK, IWORK,
167  $ PIVMIN, SPDIAM, INFO )
168 *
169 * -- LAPACK auxiliary routine --
170 * -- LAPACK is a software package provided by Univ. of Tennessee, --
171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172 *
173 * .. Scalar Arguments ..
174  INTEGER IFIRST, ILAST, INFO, N, OFFSET
175  REAL PIVMIN, RTOL, SPDIAM
176 * ..
177 * .. Array Arguments ..
178  INTEGER IWORK( * )
179  REAL D( * ), E2( * ), W( * ),
180  $ werr( * ), work( * )
181 * ..
182 *
183 * =====================================================================
184 *
185 * .. Parameters ..
186  REAL ZERO, ONE, TWO, HALF
187  PARAMETER ( ZERO = 0.0e0, one = 1.0e0, two = 2.0e0,
188  $ half = 0.5e0 )
189  INTEGER MAXITR
190 * ..
191 * .. Local Scalars ..
192  INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
193  $ OLNINT, P, PREV, SAVI1
194  REAL DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
195 *
196 * ..
197 * .. Intrinsic Functions ..
198  INTRINSIC abs, max
199 * ..
200 * .. Executable Statements ..
201 *
202  info = 0
203 *
204 * Quick return if possible
205 *
206  IF( n.LE.0 ) THEN
207  RETURN
208  END IF
209 *
210  maxitr = int( ( log( spdiam+pivmin )-log( pivmin ) ) /
211  $ log( two ) ) + 2
212 *
213 * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
214 * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
215 * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
216 * for an unconverged interval is set to the index of the next unconverged
217 * interval, and is -1 or 0 for a converged interval. Thus a linked
218 * list of unconverged intervals is set up.
219 *
220 
221  i1 = ifirst
222  i2 = ilast
223 * The number of unconverged intervals
224  nint = 0
225 * The last unconverged interval found
226  prev = 0
227  DO 75 i = i1, i2
228  k = 2*i
229  ii = i - offset
230  left = w( ii ) - werr( ii )
231  mid = w(ii)
232  right = w( ii ) + werr( ii )
233  width = right - mid
234  tmp = max( abs( left ), abs( right ) )
235 
236 * The following test prevents the test of converged intervals
237  IF( width.LT.rtol*tmp ) THEN
238 * This interval has already converged and does not need refinement.
239 * (Note that the gaps might change through refining the
240 * eigenvalues, however, they can only get bigger.)
241 * Remove it from the list.
242  iwork( k-1 ) = -1
243 * Make sure that I1 always points to the first unconverged interval
244  IF((i.EQ.i1).AND.(i.LT.i2)) i1 = i + 1
245  IF((prev.GE.i1).AND.(i.LE.i2)) iwork( 2*prev-1 ) = i + 1
246  ELSE
247 * unconverged interval found
248  prev = i
249 * Make sure that [LEFT,RIGHT] contains the desired eigenvalue
250 *
251 * Do while( CNT(LEFT).GT.I-1 )
252 *
253  fac = one
254  20 CONTINUE
255  cnt = 0
256  s = left
257  dplus = d( 1 ) - s
258  IF( dplus.LT.zero ) cnt = cnt + 1
259  DO 30 j = 2, n
260  dplus = d( j ) - s - e2( j-1 )/dplus
261  IF( dplus.LT.zero ) cnt = cnt + 1
262  30 CONTINUE
263  IF( cnt.GT.i-1 ) THEN
264  left = left - werr( ii )*fac
265  fac = two*fac
266  GO TO 20
267  END IF
268 *
269 * Do while( CNT(RIGHT).LT.I )
270 *
271  fac = one
272  50 CONTINUE
273  cnt = 0
274  s = right
275  dplus = d( 1 ) - s
276  IF( dplus.LT.zero ) cnt = cnt + 1
277  DO 60 j = 2, n
278  dplus = d( j ) - s - e2( j-1 )/dplus
279  IF( dplus.LT.zero ) cnt = cnt + 1
280  60 CONTINUE
281  IF( cnt.LT.i ) THEN
282  right = right + werr( ii )*fac
283  fac = two*fac
284  GO TO 50
285  END IF
286  nint = nint + 1
287  iwork( k-1 ) = i + 1
288  iwork( k ) = cnt
289  END IF
290  work( k-1 ) = left
291  work( k ) = right
292  75 CONTINUE
293 
294 
295  savi1 = i1
296 *
297 * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
298 * and while (ITER.LT.MAXITR)
299 *
300  iter = 0
301  80 CONTINUE
302  prev = i1 - 1
303  i = i1
304  olnint = nint
305 
306  DO 100 p = 1, olnint
307  k = 2*i
308  ii = i - offset
309  next = iwork( k-1 )
310  left = work( k-1 )
311  right = work( k )
312  mid = half*( left + right )
313 
314 * semiwidth of interval
315  width = right - mid
316  tmp = max( abs( left ), abs( right ) )
317 
318  IF( ( width.LT.rtol*tmp ) .OR.
319  $ (iter.EQ.maxitr) )THEN
320 * reduce number of unconverged intervals
321  nint = nint - 1
322 * Mark interval as converged.
323  iwork( k-1 ) = 0
324  IF( i1.EQ.i ) THEN
325  i1 = next
326  ELSE
327 * Prev holds the last unconverged interval previously examined
328  IF(prev.GE.i1) iwork( 2*prev-1 ) = next
329  END IF
330  i = next
331  GO TO 100
332  END IF
333  prev = i
334 *
335 * Perform one bisection step
336 *
337  cnt = 0
338  s = mid
339  dplus = d( 1 ) - s
340  IF( dplus.LT.zero ) cnt = cnt + 1
341  DO 90 j = 2, n
342  dplus = d( j ) - s - e2( j-1 )/dplus
343  IF( dplus.LT.zero ) cnt = cnt + 1
344  90 CONTINUE
345  IF( cnt.LE.i-1 ) THEN
346  work( k-1 ) = mid
347  ELSE
348  work( k ) = mid
349  END IF
350  i = next
351 
352  100 CONTINUE
353  iter = iter + 1
354 * do another loop if there are still unconverged intervals
355 * However, in the last iteration, all intervals are accepted
356 * since this is the best we can do.
357  IF( ( nint.GT.0 ).AND.(iter.LE.maxitr) ) GO TO 80
358 *
359 *
360 * At this point, all the intervals have converged
361  DO 110 i = savi1, ilast
362  k = 2*i
363  ii = i - offset
364 * All intervals marked by '0' have been refined.
365  IF( iwork( k-1 ).EQ.0 ) THEN
366  w( ii ) = half*( work( k-1 )+work( k ) )
367  werr( ii ) = work( k ) - w( ii )
368  END IF
369  110 CONTINUE
370 *
371 
372  RETURN
373 *
374 * End of SLARRJ
375 *
376  END
subroutine slarrj(N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO)
SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
Definition: slarrj.f:168