LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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.

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: