LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.
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 dlarrj.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 DOUBLE PRECISION PIVMIN, RTOL, SPDIAM
176* ..
177* .. Array Arguments ..
178 INTEGER IWORK( * )
179 DOUBLE PRECISION D( * ), E2( * ), W( * ),
180 $ WERR( * ), WORK( * )
181* ..
182*
183* =====================================================================
184*
185* .. Parameters ..
186 DOUBLE PRECISION ZERO, ONE, TWO, HALF
187 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
188 $ half = 0.5d0 )
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 DOUBLE PRECISION 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 DLARRJ
375*
Here is the caller graph for this function: