001:       SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
002:      $                   W, WGAP, WERR,
003:      $                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
004:      $                   DPLUS, LPLUS, WORK, INFO )
005: *
006: *  -- LAPACK auxiliary routine (version 3.2) --
007: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
008: *     November 2006
009: **
010: *     .. Scalar Arguments ..
011:       INTEGER            CLSTRT, CLEND, INFO, N
012:       DOUBLE PRECISION   CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   D( * ), DPLUS( * ), L( * ), LD( * ),
016:      $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  Given the initial representation L D L^T and its cluster of close
023: *  eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
024: *  W( CLEND ), DLARRF finds a new relatively robust representation
025: *  L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
026: *  eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
027: *
028: *  Arguments
029: *  =========
030: *
031: *  N       (input) INTEGER
032: *          The order of the matrix (subblock, if the matrix splitted).
033: *
034: *  D       (input) DOUBLE PRECISION array, dimension (N)
035: *          The N diagonal elements of the diagonal matrix D.
036: *
037: *  L       (input) DOUBLE PRECISION array, dimension (N-1)
038: *          The (N-1) subdiagonal elements of the unit bidiagonal
039: *          matrix L.
040: *
041: *  LD      (input) DOUBLE PRECISION array, dimension (N-1)
042: *          The (N-1) elements L(i)*D(i).
043: *
044: *  CLSTRT  (input) INTEGER
045: *          The index of the first eigenvalue in the cluster.
046: *
047: *  CLEND   (input) INTEGER
048: *          The index of the last eigenvalue in the cluster.
049: *
050: *  W       (input) DOUBLE PRECISION array, dimension >=  (CLEND-CLSTRT+1)
051: *          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
052: *          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
053: *          close eigenalues.
054: *
055: *  WGAP    (input/output) DOUBLE PRECISION array, dimension >=  (CLEND-CLSTRT+1)
056: *          The separation from the right neighbor eigenvalue in W.
057: *
058: *  WERR    (input) DOUBLE PRECISION array, dimension >=  (CLEND-CLSTRT+1)
059: *          WERR contain the semiwidth of the uncertainty
060: *          interval of the corresponding eigenvalue APPROXIMATION in W
061: *
062: *  SPDIAM (input) estimate of the spectral diameter obtained from the
063: *          Gerschgorin intervals
064: *
065: *  CLGAPL, CLGAPR (input) absolute gap on each end of the cluster.
066: *          Set by the calling routine to protect against shifts too close
067: *          to eigenvalues outside the cluster.
068: *
069: *  PIVMIN  (input) DOUBLE PRECISION
070: *          The minimum pivot allowed in the Sturm sequence.
071: *
072: *  SIGMA   (output) DOUBLE PRECISION
073: *          The shift used to form L(+) D(+) L(+)^T.
074: *
075: *  DPLUS   (output) DOUBLE PRECISION array, dimension (N)
076: *          The N diagonal elements of the diagonal matrix D(+).
077: *
078: *  LPLUS   (output) DOUBLE PRECISION array, dimension (N-1)
079: *          The first (N-1) elements of LPLUS contain the subdiagonal
080: *          elements of the unit bidiagonal matrix L(+).
081: *
082: *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
083: *          Workspace.
084: *
085: *  Further Details
086: *  ===============
087: *
088: *  Based on contributions by
089: *     Beresford Parlett, University of California, Berkeley, USA
090: *     Jim Demmel, University of California, Berkeley, USA
091: *     Inderjit Dhillon, University of Texas, Austin, USA
092: *     Osni Marques, LBNL/NERSC, USA
093: *     Christof Voemel, University of California, Berkeley, USA
094: *
095: *  =====================================================================
096: *
097: *     .. Parameters ..
098:       DOUBLE PRECISION   FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
099:      $                   ZERO
100:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
101:      $                     FOUR = 4.0D0, QUART = 0.25D0,
102:      $                     MAXGROWTH1 = 8.D0,
103:      $                     MAXGROWTH2 = 8.D0 )
104: *     ..
105: *     .. Local Scalars ..
106:       LOGICAL   DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
107:       INTEGER            I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
108:       PARAMETER          ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
109:       DOUBLE PRECISION   AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
110:      $                   FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
111:      $                   MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
112:      $                   RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
113: *     ..
114: *     .. External Functions ..
115:       LOGICAL DISNAN
116:       DOUBLE PRECISION   DLAMCH
117:       EXTERNAL           DISNAN, DLAMCH
118: *     ..
119: *     .. External Subroutines ..
120:       EXTERNAL           DCOPY
121: *     ..
122: *     .. Intrinsic Functions ..
123:       INTRINSIC          ABS
124: *     ..
125: *     .. Executable Statements ..
126: *
127:       INFO = 0
128:       FACT = DBLE(2**KTRYMAX)
129:       EPS = DLAMCH( 'Precision' )
130:       SHIFT = 0
131:       FORCER = .FALSE.
132: 
133: 
134: *     Note that we cannot guarantee that for any of the shifts tried,
135: *     the factorization has a small or even moderate element growth.
136: *     There could be Ritz values at both ends of the cluster and despite
137: *     backing off, there are examples where all factorizations tried
138: *     (in IEEE mode, allowing zero pivots & infinities) have INFINITE
139: *     element growth.
140: *     For this reason, we should use PIVMIN in this subroutine so that at
141: *     least the L D L^T factorization exists. It can be checked afterwards
142: *     whether the element growth caused bad residuals/orthogonality.
143: 
144: *     Decide whether the code should accept the best among all
145: *     representations despite large element growth or signal INFO=1
146:       NOFAIL = .TRUE.
147: *
148: 
149: *     Compute the average gap length of the cluster
150:       CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
151:       AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
152:       MINGAP = MIN(CLGAPL, CLGAPR)
153: *     Initial values for shifts to both ends of cluster
154:       LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
155:       RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
156: 
157: *     Use a small fudge to make sure that we really shift to the outside
158:       LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
159:       RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
160: 
161: *     Compute upper bounds for how much to back off the initial shifts
162:       LDMAX = QUART * MINGAP + TWO * PIVMIN
163:       RDMAX = QUART * MINGAP + TWO * PIVMIN
164: 
165:       LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
166:       RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
167: *
168: *     Initialize the record of the best representation found
169: *
170:       S = DLAMCH( 'S' )
171:       SMLGROWTH = ONE / S
172:       FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
173:       FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
174:       BESTSHIFT = LSIGMA
175: *
176: *     while (KTRY <= KTRYMAX)
177:       KTRY = 0
178:       GROWTHBOUND = MAXGROWTH1*SPDIAM
179: 
180:  5    CONTINUE
181:       SAWNAN1 = .FALSE.
182:       SAWNAN2 = .FALSE.
183: *     Ensure that we do not back off too much of the initial shifts
184:       LDELTA = MIN(LDMAX,LDELTA)
185:       RDELTA = MIN(RDMAX,RDELTA)
186: 
187: *     Compute the element growth when shifting to both ends of the cluster
188: *     accept the shift if there is no element growth at one of the two ends
189: 
190: *     Left end
191:       S = -LSIGMA
192:       DPLUS( 1 ) = D( 1 ) + S
193:       IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
194:          DPLUS(1) = -PIVMIN
195: *        Need to set SAWNAN1 because refined RRR test should not be used
196: *        in this case
197:          SAWNAN1 = .TRUE.
198:       ENDIF
199:       MAX1 = ABS( DPLUS( 1 ) )
200:       DO 6 I = 1, N - 1
201:          LPLUS( I ) = LD( I ) / DPLUS( I )
202:          S = S*LPLUS( I )*L( I ) - LSIGMA
203:          DPLUS( I+1 ) = D( I+1 ) + S
204:          IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
205:             DPLUS(I+1) = -PIVMIN
206: *           Need to set SAWNAN1 because refined RRR test should not be used
207: *           in this case
208:             SAWNAN1 = .TRUE.
209:          ENDIF
210:          MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
211:  6    CONTINUE
212:       SAWNAN1 = SAWNAN1 .OR.  DISNAN( MAX1 )
213: 
214:       IF( FORCER .OR.
215:      $   (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
216:          SIGMA = LSIGMA
217:          SHIFT = SLEFT
218:          GOTO 100
219:       ENDIF
220: 
221: *     Right end
222:       S = -RSIGMA
223:       WORK( 1 ) = D( 1 ) + S
224:       IF(ABS(WORK(1)).LT.PIVMIN) THEN
225:          WORK(1) = -PIVMIN
226: *        Need to set SAWNAN2 because refined RRR test should not be used
227: *        in this case
228:          SAWNAN2 = .TRUE.
229:       ENDIF
230:       MAX2 = ABS( WORK( 1 ) )
231:       DO 7 I = 1, N - 1
232:          WORK( N+I ) = LD( I ) / WORK( I )
233:          S = S*WORK( N+I )*L( I ) - RSIGMA
234:          WORK( I+1 ) = D( I+1 ) + S
235:          IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
236:             WORK(I+1) = -PIVMIN
237: *           Need to set SAWNAN2 because refined RRR test should not be used
238: *           in this case
239:             SAWNAN2 = .TRUE.
240:          ENDIF
241:          MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
242:  7    CONTINUE
243:       SAWNAN2 = SAWNAN2 .OR.  DISNAN( MAX2 )
244: 
245:       IF( FORCER .OR.
246:      $   (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
247:          SIGMA = RSIGMA
248:          SHIFT = SRIGHT
249:          GOTO 100
250:       ENDIF
251: *     If we are at this point, both shifts led to too much element growth
252: 
253: *     Record the better of the two shifts (provided it didn't lead to NaN)
254:       IF(SAWNAN1.AND.SAWNAN2) THEN
255: *        both MAX1 and MAX2 are NaN
256:          GOTO 50
257:       ELSE
258:          IF( .NOT.SAWNAN1 ) THEN
259:             INDX = 1
260:             IF(MAX1.LE.SMLGROWTH) THEN
261:                SMLGROWTH = MAX1
262:                BESTSHIFT = LSIGMA
263:             ENDIF
264:          ENDIF
265:          IF( .NOT.SAWNAN2 ) THEN
266:             IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
267:             IF(MAX2.LE.SMLGROWTH) THEN
268:                SMLGROWTH = MAX2
269:                BESTSHIFT = RSIGMA
270:             ENDIF
271:          ENDIF
272:       ENDIF
273: 
274: *     If we are here, both the left and the right shift led to
275: *     element growth. If the element growth is moderate, then
276: *     we may still accept the representation, if it passes a
277: *     refined test for RRR. This test supposes that no NaN occurred.
278: *     Moreover, we use the refined RRR test only for isolated clusters.
279:       IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
280:      $   (MIN(MAX1,MAX2).LT.FAIL2)
281:      $  .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
282:          DORRR1 = .TRUE.
283:       ELSE
284:          DORRR1 = .FALSE.
285:       ENDIF
286:       TRYRRR1 = .TRUE.
287:       IF( TRYRRR1 .AND. DORRR1 ) THEN
288:       IF(INDX.EQ.1) THEN
289:          TMP = ABS( DPLUS( N ) )
290:          ZNM2 = ONE
291:          PROD = ONE
292:          OLDP = ONE
293:          DO 15 I = N-1, 1, -1
294:             IF( PROD .LE. EPS ) THEN
295:                PROD =
296:      $         ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
297:             ELSE
298:                PROD = PROD*ABS(WORK(N+I))
299:             END IF
300:             OLDP = PROD
301:             ZNM2 = ZNM2 + PROD**2
302:             TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
303:  15      CONTINUE
304:          RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
305:          IF (RRR1.LE.MAXGROWTH2) THEN
306:             SIGMA = LSIGMA
307:             SHIFT = SLEFT
308:             GOTO 100
309:          ENDIF
310:       ELSE IF(INDX.EQ.2) THEN
311:          TMP = ABS( WORK( N ) )
312:          ZNM2 = ONE
313:          PROD = ONE
314:          OLDP = ONE
315:          DO 16 I = N-1, 1, -1
316:             IF( PROD .LE. EPS ) THEN
317:                PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
318:             ELSE
319:                PROD = PROD*ABS(LPLUS(I))
320:             END IF
321:             OLDP = PROD
322:             ZNM2 = ZNM2 + PROD**2
323:             TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
324:  16      CONTINUE
325:          RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
326:          IF (RRR2.LE.MAXGROWTH2) THEN
327:             SIGMA = RSIGMA
328:             SHIFT = SRIGHT
329:             GOTO 100
330:          ENDIF
331:       END IF
332:       ENDIF
333: 
334:  50   CONTINUE
335: 
336:       IF (KTRY.LT.KTRYMAX) THEN
337: *        If we are here, both shifts failed also the RRR test.
338: *        Back off to the outside
339:          LSIGMA = MAX( LSIGMA - LDELTA,
340:      $     LSIGMA - LDMAX)
341:          RSIGMA = MIN( RSIGMA + RDELTA,
342:      $     RSIGMA + RDMAX )
343:          LDELTA = TWO * LDELTA
344:          RDELTA = TWO * RDELTA
345:          KTRY = KTRY + 1
346:          GOTO 5
347:       ELSE
348: *        None of the representations investigated satisfied our
349: *        criteria. Take the best one we found.
350:          IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
351:             LSIGMA = BESTSHIFT
352:             RSIGMA = BESTSHIFT
353:             FORCER = .TRUE.
354:             GOTO 5
355:          ELSE
356:             INFO = 1
357:             RETURN
358:          ENDIF
359:       END IF
360: 
361:  100  CONTINUE
362:       IF (SHIFT.EQ.SLEFT) THEN
363:       ELSEIF (SHIFT.EQ.SRIGHT) THEN
364: *        store new L and D back into DPLUS, LPLUS
365:          CALL DCOPY( N, WORK, 1, DPLUS, 1 )
366:          CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
367:       ENDIF
368: 
369:       RETURN
370: *
371: *     End of DLARRF
372: *
373:       END
374: