LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlarrf.f
Go to the documentation of this file.
1*> \brief \b DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is relatively isolated.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLARRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
22* W, WGAP, WERR,
23* SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
24* DPLUS, LPLUS, WORK, INFO )
25*
26* .. Scalar Arguments ..
27* INTEGER CLSTRT, CLEND, INFO, N
28* DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
29* ..
30* .. Array Arguments ..
31* DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
32* $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> Given the initial representation L D L^T and its cluster of close
42*> eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
43*> W( CLEND ), DLARRF finds a new relatively robust representation
44*> L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
45*> eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] N
52*> \verbatim
53*> N is INTEGER
54*> The order of the matrix (subblock, if the matrix split).
55*> \endverbatim
56*>
57*> \param[in] D
58*> \verbatim
59*> D is DOUBLE PRECISION array, dimension (N)
60*> The N diagonal elements of the diagonal matrix D.
61*> \endverbatim
62*>
63*> \param[in] L
64*> \verbatim
65*> L is DOUBLE PRECISION array, dimension (N-1)
66*> The (N-1) subdiagonal elements of the unit bidiagonal
67*> matrix L.
68*> \endverbatim
69*>
70*> \param[in] LD
71*> \verbatim
72*> LD is DOUBLE PRECISION array, dimension (N-1)
73*> The (N-1) elements L(i)*D(i).
74*> \endverbatim
75*>
76*> \param[in] CLSTRT
77*> \verbatim
78*> CLSTRT is INTEGER
79*> The index of the first eigenvalue in the cluster.
80*> \endverbatim
81*>
82*> \param[in] CLEND
83*> \verbatim
84*> CLEND is INTEGER
85*> The index of the last eigenvalue in the cluster.
86*> \endverbatim
87*>
88*> \param[in] W
89*> \verbatim
90*> W is DOUBLE PRECISION array, dimension
91*> dimension is >= (CLEND-CLSTRT+1)
92*> The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
93*> W( CLSTRT ) through W( CLEND ) form the cluster of relatively
94*> close eigenalues.
95*> \endverbatim
96*>
97*> \param[in,out] WGAP
98*> \verbatim
99*> WGAP is DOUBLE PRECISION array, dimension
100*> dimension is >= (CLEND-CLSTRT+1)
101*> The separation from the right neighbor eigenvalue in W.
102*> \endverbatim
103*>
104*> \param[in] WERR
105*> \verbatim
106*> WERR is DOUBLE PRECISION array, dimension
107*> dimension is >= (CLEND-CLSTRT+1)
108*> WERR contain the semiwidth of the uncertainty
109*> interval of the corresponding eigenvalue APPROXIMATION in W
110*> \endverbatim
111*>
112*> \param[in] SPDIAM
113*> \verbatim
114*> SPDIAM is DOUBLE PRECISION
115*> estimate of the spectral diameter obtained from the
116*> Gerschgorin intervals
117*> \endverbatim
118*>
119*> \param[in] CLGAPL
120*> \verbatim
121*> CLGAPL is DOUBLE PRECISION
122*> \endverbatim
123*>
124*> \param[in] CLGAPR
125*> \verbatim
126*> CLGAPR is DOUBLE PRECISION
127*> absolute gap on each end of the cluster.
128*> Set by the calling routine to protect against shifts too close
129*> to eigenvalues outside the cluster.
130*> \endverbatim
131*>
132*> \param[in] PIVMIN
133*> \verbatim
134*> PIVMIN is DOUBLE PRECISION
135*> The minimum pivot allowed in the Sturm sequence.
136*> \endverbatim
137*>
138*> \param[out] SIGMA
139*> \verbatim
140*> SIGMA is DOUBLE PRECISION
141*> The shift used to form L(+) D(+) L(+)^T.
142*> \endverbatim
143*>
144*> \param[out] DPLUS
145*> \verbatim
146*> DPLUS is DOUBLE PRECISION array, dimension (N)
147*> The N diagonal elements of the diagonal matrix D(+).
148*> \endverbatim
149*>
150*> \param[out] LPLUS
151*> \verbatim
152*> LPLUS is DOUBLE PRECISION array, dimension (N-1)
153*> The first (N-1) elements of LPLUS contain the subdiagonal
154*> elements of the unit bidiagonal matrix L(+).
155*> \endverbatim
156*>
157*> \param[out] WORK
158*> \verbatim
159*> WORK is DOUBLE PRECISION array, dimension (2*N)
160*> Workspace.
161*> \endverbatim
162*>
163*> \param[out] INFO
164*> \verbatim
165*> INFO is INTEGER
166*> Signals processing OK (=0) or failure (=1)
167*> \endverbatim
168*
169* Authors:
170* ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup larrf
178*
179*> \par Contributors:
180* ==================
181*>
182*> Beresford Parlett, University of California, Berkeley, USA \n
183*> Jim Demmel, University of California, Berkeley, USA \n
184*> Inderjit Dhillon, University of Texas, Austin, USA \n
185*> Osni Marques, LBNL/NERSC, USA \n
186*> Christof Voemel, University of California, Berkeley, USA
187*
188* =====================================================================
189 SUBROUTINE dlarrf( N, D, L, LD, CLSTRT, CLEND,
190 $ W, WGAP, WERR,
191 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
192 $ DPLUS, LPLUS, WORK, INFO )
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 DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
201* ..
202* .. Array Arguments ..
203 DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
204 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO
211 PARAMETER ( ONE = 1.0d0, two = 2.0d0, four = 4.0d0,
212 $ quart = 0.25d0,
213 $ maxgrowth1 = 8.d0,
214 $ maxgrowth2 = 8.d0 )
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 DOUBLE PRECISION 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 DISNAN
227 DOUBLE PRECISION DLAMCH
228 EXTERNAL DISNAN, DLAMCH
229* ..
230* .. External Subroutines ..
231 EXTERNAL dcopy
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 = dble(2**ktrymax)
247 eps = dlamch( '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 / dble(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)* four * eps
278 rsigma = rsigma + abs(rsigma)* four * 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 = dlamch( 'S' )
290 smlgrowth = one / s
291 fail = dble(n-1)*mingap/(spdiam*eps)
292 fail2 = dble(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. disnan( 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. disnan( 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/dble(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 dcopy( n, work, 1, dplus, 1 )
485 CALL dcopy( n-1, work(n+1), 1, lplus, 1 )
486 ENDIF
487
488 RETURN
489*
490* End of DLARRF
491*
492 END
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlarrf(n, d, l, ld, clstrt, clend, w, wgap, werr, spdiam, clgapl, clgapr, pivmin, sigma, dplus, lplus, work, info)
DLARRF finds a new relatively robust representation such that at least one of the eigenvalues is rela...
Definition dlarrf.f:193