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*> \date June 2016
178*
179*> \ingroup OTHERauxiliary
180*
181*> \par Contributors:
182*  ==================
183*>
184*> Beresford Parlett, University of California, Berkeley, USA \n
185*> Jim Demmel, University of California, Berkeley, USA \n
186*> Inderjit Dhillon, University of Texas, Austin, USA \n
187*> Osni Marques, LBNL/NERSC, USA \n
188*> Christof Voemel, University of California, Berkeley, USA
189*
190*  =====================================================================
191      SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
192     $                   W, WGAP, WERR,
193     $                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
194     $                   DPLUS, LPLUS, WORK, INFO )
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*
495      END
496