LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2016
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 195 of file dlarrf.f.

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