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