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

◆ dlarrf()

subroutine dlarrf ( integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) l,
double precision, dimension( * ) ld,
integer clstrt,
integer clend,
double precision, dimension( * ) w,
double precision, dimension( * ) wgap,
double precision, dimension( * ) werr,
double precision spdiam,
double precision clgapl,
double precision clgapr,
double precision pivmin,
double precision sigma,
double precision, dimension( * ) dplus,
double precision, dimension( * ) lplus,
double precision, dimension( * ) work,
integer info )

DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.

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

Purpose:
!>
!> Given the initial representation L D L^T and its cluster of close
!> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
!> W( CLEND ), DLARRF finds a new relatively robust representation
!> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
!> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix (subblock, if the matrix split).
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The N diagonal elements of the diagonal matrix D.
!>
[in]L
!>          L is DOUBLE PRECISION array, dimension (N-1)
!>          The (N-1) subdiagonal elements of the unit bidiagonal
!>          matrix L.
!>
[in]LD
!>          LD is DOUBLE PRECISION array, dimension (N-1)
!>          The (N-1) elements L(i)*D(i).
!>
[in]CLSTRT
!>          CLSTRT is INTEGER
!>          The index of the first eigenvalue in the cluster.
!>
[in]CLEND
!>          CLEND is INTEGER
!>          The index of the last eigenvalue in the cluster.
!>
[in]W
!>          W is DOUBLE PRECISION array, dimension
!>          dimension is >=  (CLEND-CLSTRT+1)
!>          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
!>          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
!>          close eigenalues.
!>
[in,out]WGAP
!>          WGAP is DOUBLE PRECISION array, dimension
!>          dimension is >=  (CLEND-CLSTRT+1)
!>          The separation from the right neighbor eigenvalue in W.
!>
[in]WERR
!>          WERR is DOUBLE PRECISION array, dimension
!>          dimension is  >=  (CLEND-CLSTRT+1)
!>          WERR contain the semiwidth of the uncertainty
!>          interval of the corresponding eigenvalue APPROXIMATION in W
!>
[in]SPDIAM
!>          SPDIAM is DOUBLE PRECISION
!>          estimate of the spectral diameter obtained from the
!>          Gerschgorin intervals
!>
[in]CLGAPL
!>          CLGAPL is DOUBLE PRECISION
!>
[in]CLGAPR
!>          CLGAPR is DOUBLE PRECISION
!>          absolute gap on each end of the cluster.
!>          Set by the calling routine to protect against shifts too close
!>          to eigenvalues outside the cluster.
!>
[in]PIVMIN
!>          PIVMIN is DOUBLE PRECISION
!>          The minimum pivot allowed in the Sturm sequence.
!>
[out]SIGMA
!>          SIGMA is DOUBLE PRECISION
!>          The shift used to form L(+) D(+) L(+)^T.
!>
[out]DPLUS
!>          DPLUS is DOUBLE PRECISION array, dimension (N)
!>          The N diagonal elements of the diagonal matrix D(+).
!>
[out]LPLUS
!>          LPLUS is DOUBLE PRECISION array, dimension (N-1)
!>          The first (N-1) elements of LPLUS contain the subdiagonal
!>          elements of the unit bidiagonal matrix L(+).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*N)
!>          Workspace.
!>
[out]INFO
!>          INFO is INTEGER
!>          Signals processing OK (=0) or failure (=1)
!>
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 187 of file dlarrf.f.

191*
192* -- LAPACK auxiliary routine --
193* -- LAPACK is a software package provided by Univ. of Tennessee, --
194* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*
196* .. Scalar Arguments ..
197 INTEGER CLSTRT, CLEND, INFO, N
198 DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
199* ..
200* .. Array Arguments ..
201 DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
202 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
203* ..
204*
205* =====================================================================
206*
207* .. Parameters ..
208 DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO
209 parameter( one = 1.0d0, two = 2.0d0, four = 4.0d0,
210 $ quart = 0.25d0,
211 $ maxgrowth1 = 8.d0,
212 $ maxgrowth2 = 8.d0 )
213* ..
214* .. Local Scalars ..
215 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
216 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
217 parameter( ktrymax = 1, sleft = 1, sright = 2 )
218 DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
219 $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
220 $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
221 $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
222* ..
223* .. External Functions ..
224 LOGICAL DISNAN
225 DOUBLE PRECISION DLAMCH
226 EXTERNAL disnan, dlamch
227* ..
228* .. External Subroutines ..
229 EXTERNAL dcopy
230* ..
231* .. Intrinsic Functions ..
232 INTRINSIC abs
233* ..
234* .. Executable Statements ..
235*
236 info = 0
237*
238* Quick return if possible
239*
240 IF( n.LE.0 ) THEN
241 RETURN
242 END IF
243*
244 fact = dble(2**ktrymax)
245 eps = dlamch( 'Precision' )
246 shift = 0
247 forcer = .false.
248
249
250* Note that we cannot guarantee that for any of the shifts tried,
251* the factorization has a small or even moderate element growth.
252* There could be Ritz values at both ends of the cluster and despite
253* backing off, there are examples where all factorizations tried
254* (in IEEE mode, allowing zero pivots & infinities) have INFINITE
255* element growth.
256* For this reason, we should use PIVMIN in this subroutine so that at
257* least the L D L^T factorization exists. It can be checked afterwards
258* whether the element growth caused bad residuals/orthogonality.
259
260* Decide whether the code should accept the best among all
261* representations despite large element growth or signal INFO=1
262* Setting NOFAIL to .FALSE. for quick fix for bug 113
263 nofail = .false.
264*
265
266* Compute the average gap length of the cluster
267 clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
268 avgap = clwdth / dble(clend-clstrt)
269 mingap = min(clgapl, clgapr)
270* Initial values for shifts to both ends of cluster
271 lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
272 rsigma = max(w( clstrt ),w( clend )) + werr( clend )
273
274* Use a small fudge to make sure that we really shift to the outside
275 lsigma = lsigma - abs(lsigma)* four * eps
276 rsigma = rsigma + abs(rsigma)* four * eps
277
278* Compute upper bounds for how much to back off the initial shifts
279 ldmax = quart * mingap + two * pivmin
280 rdmax = quart * mingap + two * pivmin
281
282 ldelta = max(avgap,wgap( clstrt ))/fact
283 rdelta = max(avgap,wgap( clend-1 ))/fact
284*
285* Initialize the record of the best representation found
286*
287 s = dlamch( 'S' )
288 smlgrowth = one / s
289 fail = dble(n-1)*mingap/(spdiam*eps)
290 fail2 = dble(n-1)*mingap/(spdiam*sqrt(eps))
291 bestshift = lsigma
292*
293* while (KTRY <= KTRYMAX)
294 ktry = 0
295 growthbound = maxgrowth1*spdiam
296
297 5 CONTINUE
298 sawnan1 = .false.
299 sawnan2 = .false.
300* Ensure that we do not back off too much of the initial shifts
301 ldelta = min(ldmax,ldelta)
302 rdelta = min(rdmax,rdelta)
303
304* Compute the element growth when shifting to both ends of the cluster
305* accept the shift if there is no element growth at one of the two ends
306
307* Left end
308 s = -lsigma
309 dplus( 1 ) = d( 1 ) + s
310 IF(abs(dplus(1)).LT.pivmin) THEN
311 dplus(1) = -pivmin
312* Need to set SAWNAN1 because refined RRR test should not be used
313* in this case
314 sawnan1 = .true.
315 ENDIF
316 max1 = abs( dplus( 1 ) )
317 DO 6 i = 1, n - 1
318 lplus( i ) = ld( i ) / dplus( i )
319 s = s*lplus( i )*l( i ) - lsigma
320 dplus( i+1 ) = d( i+1 ) + s
321 IF(abs(dplus(i+1)).LT.pivmin) THEN
322 dplus(i+1) = -pivmin
323* Need to set SAWNAN1 because refined RRR test should not be used
324* in this case
325 sawnan1 = .true.
326 ENDIF
327 max1 = max( max1,abs(dplus(i+1)) )
328 6 CONTINUE
329 sawnan1 = sawnan1 .OR. disnan( max1 )
330
331 IF( forcer .OR.
332 $ (max1.LE.growthbound .AND. .NOT.sawnan1 ) ) THEN
333 sigma = lsigma
334 shift = sleft
335 GOTO 100
336 ENDIF
337
338* Right end
339 s = -rsigma
340 work( 1 ) = d( 1 ) + s
341 IF(abs(work(1)).LT.pivmin) THEN
342 work(1) = -pivmin
343* Need to set SAWNAN2 because refined RRR test should not be used
344* in this case
345 sawnan2 = .true.
346 ENDIF
347 max2 = abs( work( 1 ) )
348 DO 7 i = 1, n - 1
349 work( n+i ) = ld( i ) / work( i )
350 s = s*work( n+i )*l( i ) - rsigma
351 work( i+1 ) = d( i+1 ) + s
352 IF(abs(work(i+1)).LT.pivmin) THEN
353 work(i+1) = -pivmin
354* Need to set SAWNAN2 because refined RRR test should not be used
355* in this case
356 sawnan2 = .true.
357 ENDIF
358 max2 = max( max2,abs(work(i+1)) )
359 7 CONTINUE
360 sawnan2 = sawnan2 .OR. disnan( max2 )
361
362 IF( forcer .OR.
363 $ (max2.LE.growthbound .AND. .NOT.sawnan2 ) ) THEN
364 sigma = rsigma
365 shift = sright
366 GOTO 100
367 ENDIF
368* If we are at this point, both shifts led to too much element growth
369
370* Record the better of the two shifts (provided it didn't lead to NaN)
371 IF(sawnan1.AND.sawnan2) THEN
372* both MAX1 and MAX2 are NaN
373 GOTO 50
374 ELSE
375 IF( .NOT.sawnan1 ) THEN
376 indx = 1
377 IF(max1.LE.smlgrowth) THEN
378 smlgrowth = max1
379 bestshift = lsigma
380 ENDIF
381 ENDIF
382 IF( .NOT.sawnan2 ) THEN
383 IF(sawnan1 .OR. max2.LE.max1) indx = 2
384 IF(max2.LE.smlgrowth) THEN
385 smlgrowth = max2
386 bestshift = rsigma
387 ENDIF
388 ENDIF
389 ENDIF
390
391* If we are here, both the left and the right shift led to
392* element growth. If the element growth is moderate, then
393* we may still accept the representation, if it passes a
394* refined test for RRR. This test supposes that no NaN occurred.
395* Moreover, we use the refined RRR test only for isolated clusters.
396 IF((clwdth.LT.mingap/dble(128)) .AND.
397 $ (min(max1,max2).LT.fail2)
398 $ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2)) THEN
399 dorrr1 = .true.
400 ELSE
401 dorrr1 = .false.
402 ENDIF
403 tryrrr1 = .true.
404 IF( tryrrr1 .AND. dorrr1 ) THEN
405 IF(indx.EQ.1) THEN
406 tmp = abs( dplus( n ) )
407 znm2 = one
408 prod = one
409 oldp = one
410 DO 15 i = n-1, 1, -1
411 IF( prod .LE. eps ) THEN
412 prod =
413 $ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
414 ELSE
415 prod = prod*abs(work(n+i))
416 END IF
417 oldp = prod
418 znm2 = znm2 + prod**2
419 tmp = max( tmp, abs( dplus( i ) * prod ))
420 15 CONTINUE
421 rrr1 = tmp/( spdiam * sqrt( znm2 ) )
422 IF (rrr1.LE.maxgrowth2) THEN
423 sigma = lsigma
424 shift = sleft
425 GOTO 100
426 ENDIF
427 ELSE IF(indx.EQ.2) THEN
428 tmp = abs( work( n ) )
429 znm2 = one
430 prod = one
431 oldp = one
432 DO 16 i = n-1, 1, -1
433 IF( prod .LE. eps ) THEN
434 prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
435 ELSE
436 prod = prod*abs(lplus(i))
437 END IF
438 oldp = prod
439 znm2 = znm2 + prod**2
440 tmp = max( tmp, abs( work( i ) * prod ))
441 16 CONTINUE
442 rrr2 = tmp/( spdiam * sqrt( znm2 ) )
443 IF (rrr2.LE.maxgrowth2) THEN
444 sigma = rsigma
445 shift = sright
446 GOTO 100
447 ENDIF
448 END IF
449 ENDIF
450
451 50 CONTINUE
452
453 IF (ktry.LT.ktrymax) THEN
454* If we are here, both shifts failed also the RRR test.
455* Back off to the outside
456 lsigma = max( lsigma - ldelta,
457 $ lsigma - ldmax)
458 rsigma = min( rsigma + rdelta,
459 $ rsigma + rdmax )
460 ldelta = two * ldelta
461 rdelta = two * rdelta
462 ktry = ktry + 1
463 GOTO 5
464 ELSE
465* None of the representations investigated satisfied our
466* criteria. Take the best one we found.
467 IF((smlgrowth.LT.fail).OR.nofail) THEN
468 lsigma = bestshift
469 rsigma = bestshift
470 forcer = .true.
471 GOTO 5
472 ELSE
473 info = 1
474 RETURN
475 ENDIF
476 END IF
477
478 100 CONTINUE
479 IF (shift.EQ.sleft) THEN
480 ELSEIF (shift.EQ.sright) THEN
481* store new L and D back into DPLUS, LPLUS
482 CALL dcopy( n, work, 1, dplus, 1 )
483 CALL dcopy( n-1, work(n+1), 1, lplus, 1 )
484 ENDIF
485
486 RETURN
487*
488* End of DLARRF
489*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:57
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
Here is the call graph for this function:
Here is the caller graph for this function: