LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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 split).``` [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)```
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 slarrf.f.

195 *
196 * -- LAPACK auxiliary routine (version 3.6.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  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 * Setting NOFAIL to .FALSE. for quick fix for bug 113
261  nofail = .false.
262 *
263
264 * Compute the average gap length of the cluster
265  clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
266  avgap = clwdth / REAL(clend-clstrt)
267  mingap = min(clgapl, clgapr)
268 * Initial values for shifts to both ends of cluster
269  lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
270  rsigma = max(w( clstrt ),w( clend )) + werr( clend )
271
272 * Use a small fudge to make sure that we really shift to the outside
273  lsigma = lsigma - abs(lsigma)* two * eps
274  rsigma = rsigma + abs(rsigma)* two * eps
275
276 * Compute upper bounds for how much to back off the initial shifts
277  ldmax = quart * mingap + two * pivmin
278  rdmax = quart * mingap + two * pivmin
279
280  ldelta = max(avgap,wgap( clstrt ))/fact
281  rdelta = max(avgap,wgap( clend-1 ))/fact
282 *
283 * Initialize the record of the best representation found
284 *
285  s = slamch( 'S' )
286  smlgrowth = one / s
287  fail = REAL(n-1)*mingap/(spdiam*eps)
288  fail2 = REAL(n-1)*mingap/(spdiam*sqrt(eps))
289  bestshift = lsigma
290 *
291 * while (KTRY <= KTRYMAX)
292  ktry = 0
293  growthbound = maxgrowth1*spdiam
294
295  5 CONTINUE
296  sawnan1 = .false.
297  sawnan2 = .false.
298 * Ensure that we do not back off too much of the initial shifts
299  ldelta = min(ldmax,ldelta)
300  rdelta = min(rdmax,rdelta)
301
302 * Compute the element growth when shifting to both ends of the cluster
303 * accept the shift if there is no element growth at one of the two ends
304
305 * Left end
306  s = -lsigma
307  dplus( 1 ) = d( 1 ) + s
308  IF(abs(dplus(1)).LT.pivmin) THEN
309  dplus(1) = -pivmin
310 * Need to set SAWNAN1 because refined RRR test should not be used
311 * in this case
312  sawnan1 = .true.
313  ENDIF
314  max1 = abs( dplus( 1 ) )
315  DO 6 i = 1, n - 1
316  lplus( i ) = ld( i ) / dplus( i )
317  s = s*lplus( i )*l( i ) - lsigma
318  dplus( i+1 ) = d( i+1 ) + s
319  IF(abs(dplus(i+1)).LT.pivmin) THEN
320  dplus(i+1) = -pivmin
321 * Need to set SAWNAN1 because refined RRR test should not be used
322 * in this case
323  sawnan1 = .true.
324  ENDIF
325  max1 = max( max1,abs(dplus(i+1)) )
326  6 CONTINUE
327  sawnan1 = sawnan1 .OR. sisnan( max1 )
328
329  IF( forcer .OR.
330  \$ (max1.LE.growthbound .AND. .NOT.sawnan1 ) ) THEN
331  sigma = lsigma
332  shift = sleft
333  GOTO 100
334  ENDIF
335
336 * Right end
337  s = -rsigma
338  work( 1 ) = d( 1 ) + s
339  IF(abs(work(1)).LT.pivmin) THEN
340  work(1) = -pivmin
341 * Need to set SAWNAN2 because refined RRR test should not be used
342 * in this case
343  sawnan2 = .true.
344  ENDIF
345  max2 = abs( work( 1 ) )
346  DO 7 i = 1, n - 1
347  work( n+i ) = ld( i ) / work( i )
348  s = s*work( n+i )*l( i ) - rsigma
349  work( i+1 ) = d( i+1 ) + s
350  IF(abs(work(i+1)).LT.pivmin) THEN
351  work(i+1) = -pivmin
352 * Need to set SAWNAN2 because refined RRR test should not be used
353 * in this case
354  sawnan2 = .true.
355  ENDIF
356  max2 = max( max2,abs(work(i+1)) )
357  7 CONTINUE
358  sawnan2 = sawnan2 .OR. sisnan( max2 )
359
360  IF( forcer .OR.
361  \$ (max2.LE.growthbound .AND. .NOT.sawnan2 ) ) THEN
362  sigma = rsigma
363  shift = sright
364  GOTO 100
365  ENDIF
366 * If we are at this point, both shifts led to too much element growth
367
368 * Record the better of the two shifts (provided it didn't lead to NaN)
369  IF(sawnan1.AND.sawnan2) THEN
370 * both MAX1 and MAX2 are NaN
371  GOTO 50
372  ELSE
373  IF( .NOT.sawnan1 ) THEN
374  indx = 1
375  IF(max1.LE.smlgrowth) THEN
376  smlgrowth = max1
377  bestshift = lsigma
378  ENDIF
379  ENDIF
380  IF( .NOT.sawnan2 ) THEN
381  IF(sawnan1 .OR. max2.LE.max1) indx = 2
382  IF(max2.LE.smlgrowth) THEN
383  smlgrowth = max2
384  bestshift = rsigma
385  ENDIF
386  ENDIF
387  ENDIF
388
389 * If we are here, both the left and the right shift led to
390 * element growth. If the element growth is moderate, then
391 * we may still accept the representation, if it passes a
392 * refined test for RRR. This test supposes that no NaN occurred.
393 * Moreover, we use the refined RRR test only for isolated clusters.
394  IF((clwdth.LT.mingap/REAL(128)) .AND.
395  \$ (min(max1,max2).LT.fail2)
396  \$ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2)) THEN
397  dorrr1 = .true.
398  ELSE
399  dorrr1 = .false.
400  ENDIF
401  tryrrr1 = .true.
402  IF( tryrrr1 .AND. dorrr1 ) THEN
403  IF(indx.EQ.1) THEN
404  tmp = abs( dplus( n ) )
405  znm2 = one
406  prod = one
407  oldp = one
408  DO 15 i = n-1, 1, -1
409  IF( prod .LE. eps ) THEN
410  prod =
411  \$ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
412  ELSE
413  prod = prod*abs(work(n+i))
414  END IF
415  oldp = prod
416  znm2 = znm2 + prod**2
417  tmp = max( tmp, abs( dplus( i ) * prod ))
418  15 CONTINUE
419  rrr1 = tmp/( spdiam * sqrt( znm2 ) )
420  IF (rrr1.LE.maxgrowth2) THEN
421  sigma = lsigma
422  shift = sleft
423  GOTO 100
424  ENDIF
425  ELSE IF(indx.EQ.2) THEN
426  tmp = abs( work( n ) )
427  znm2 = one
428  prod = one
429  oldp = one
430  DO 16 i = n-1, 1, -1
431  IF( prod .LE. eps ) THEN
432  prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
433  ELSE
434  prod = prod*abs(lplus(i))
435  END IF
436  oldp = prod
437  znm2 = znm2 + prod**2
438  tmp = max( tmp, abs( work( i ) * prod ))
439  16 CONTINUE
440  rrr2 = tmp/( spdiam * sqrt( znm2 ) )
441  IF (rrr2.LE.maxgrowth2) THEN
442  sigma = rsigma
443  shift = sright
444  GOTO 100
445  ENDIF
446  END IF
447  ENDIF
448
449  50 CONTINUE
450
451  IF (ktry.LT.ktrymax) THEN
452 * If we are here, both shifts failed also the RRR test.
453 * Back off to the outside
454  lsigma = max( lsigma - ldelta,
455  \$ lsigma - ldmax)
456  rsigma = min( rsigma + rdelta,
457  \$ rsigma + rdmax )
458  ldelta = two * ldelta
459  rdelta = two * rdelta
460  ktry = ktry + 1
461  GOTO 5
462  ELSE
463 * None of the representations investigated satisfied our
464 * criteria. Take the best one we found.
465  IF((smlgrowth.LT.fail).OR.nofail) THEN
466  lsigma = bestshift
467  rsigma = bestshift
468  forcer = .true.
469  GOTO 5
470  ELSE
471  info = 1
472  RETURN
473  ENDIF
474  END IF
475
476  100 CONTINUE
477  IF (shift.EQ.sleft) THEN
478  ELSEIF (shift.EQ.sright) THEN
479 * store new L and D back into DPLUS, LPLUS
480  CALL scopy( n, work, 1, dplus, 1 )
481  CALL scopy( n-1, work(n+1), 1, lplus, 1 )
482  ENDIF
483
484  RETURN
485 *
486 * End of SLARRF
487 *
logical function sisnan(SIN)
SISNAN tests input for NaN.
Definition: sisnan.f:61
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: