LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
slarrf.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine slarrf (N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO)
 SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated. More...
 

Function/Subroutine Documentation

subroutine slarrf ( integer  N,
real, dimension( * )  D,
real, dimension( * )  L,
real, dimension( * )  LD,
integer  CLSTRT,
integer  CLEND,
real, dimension( * )  W,
real, dimension( * )  WGAP,
real, dimension( * )  WERR,
real  SPDIAM,
real  CLGAPL,
real  CLGAPR,
real  PIVMIN,
real  SIGMA,
real, dimension( * )  DPLUS,
real, dimension( * )  LPLUS,
real, dimension( * )  WORK,
integer  INFO 
)

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

Download SLARRF + 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 ), SLARRF 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 splitted).
[in]D
          D is REAL array, dimension (N)
          The N diagonal elements of the diagonal matrix D.
[in]L
          L is REAL array, dimension (N-1)
          The (N-1) subdiagonal elements of the unit bidiagonal
          matrix L.
[in]LD
          LD is REAL 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 REAL 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 REAL array, dimension
          dimension is >=  (CLEND-CLSTRT+1)
          The separation from the right neighbor eigenvalue in W.
[in]WERR
          WERR is REAL 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 REAL
          estimate of the spectral diameter obtained from the
          Gerschgorin intervals
[in]CLGAPL
          CLGAPL is REAL
[in]CLGAPR
          CLGAPR is REAL
          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 REAL
          The minimum pivot allowed in the Sturm sequence.
[out]SIGMA
          SIGMA is REAL
          The shift used to form L(+) D(+) L(+)^T.
[out]DPLUS
          DPLUS is REAL array, dimension (N)
          The N diagonal elements of the diagonal matrix D(+).
[out]LPLUS
          LPLUS is REAL 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 REAL 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
September 2012
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 slarrf.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: