LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slarrf()

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)
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 189 of file slarrf.f.

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