LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dlarrj()

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

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

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

Purpose:
 Given the initial eigenvalue approximations of T, DLARRJ
 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 DOUBLE PRECISION array, dimension (N)
          The N diagonal elements of T.
[in]E2
          E2 is DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (2*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (2*N)
          Workspace.
[in]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence for T.
[in]SPDIAM
          SPDIAM is DOUBLE PRECISION
          The spectral diameter of T.
[out]INFO
          INFO is INTEGER
          Error flag.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2017
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 170 of file dlarrj.f.

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