 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ slarrj()

 subroutine slarrj ( integer N, real, dimension( * ) D, real, dimension( * ) E2, integer IFIRST, integer ILAST, real RTOL, integer OFFSET, real, dimension( * ) W, real, dimension( * ) WERR, real, dimension( * ) WORK, integer, dimension( * ) IWORK, real PIVMIN, real SPDIAM, integer INFO )

SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.

Download SLARRJ + dependencies [TGZ] [ZIP] [TXT]

Purpose:
``` Given the initial eigenvalue approximations of T, SLARRJ
does  bisection to refine the eigenvalues of T,
W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
guesses for these eigenvalues are input in W, the corresponding estimate
of the error in these guesses in WERR. During bisection, intervals
[left, right] are maintained by storing their mid-points and
semi-widths in the arrays W and WERR respectively.```
Parameters
 [in] N ``` N is INTEGER The order of the matrix.``` [in] D ``` D is REAL array, dimension (N) The N diagonal elements of T.``` [in] E2 ``` E2 is REAL array, dimension (N-1) The Squares of the (N-1) subdiagonal elements of T.``` [in] IFIRST ``` IFIRST is INTEGER The index of the first eigenvalue to be computed.``` [in] ILAST ``` ILAST is INTEGER The index of the last eigenvalue to be computed.``` [in] RTOL ``` RTOL is REAL Tolerance for the convergence of the bisection intervals. An interval [LEFT,RIGHT] has converged if RIGHT-LEFT < RTOL*MAX(|LEFT|,|RIGHT|).``` [in] OFFSET ``` OFFSET is INTEGER Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET through ILAST-OFFSET elements of these arrays are to be used.``` [in,out] W ``` W is REAL array, dimension (N) On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are estimates of the eigenvalues of L D L^T indexed IFIRST through ILAST. On output, these estimates are refined.``` [in,out] WERR ``` WERR is REAL array, dimension (N) On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are the errors in the estimates of the corresponding elements in W. On output, these errors are refined.``` [out] WORK ``` WORK is REAL array, dimension (2*N) Workspace.``` [out] IWORK ``` IWORK is INTEGER array, dimension (2*N) Workspace.``` [in] PIVMIN ``` PIVMIN is REAL The minimum pivot in the Sturm sequence for T.``` [in] SPDIAM ``` SPDIAM is REAL The spectral diameter of T.``` [out] INFO ``` INFO is INTEGER Error flag.```
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 165 of file slarrj.f.

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 *
Here is the caller graph for this function: