001:       SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST,
002:      $                   RTOL, OFFSET, W, WERR, WORK, IWORK,
003:      $                   PIVMIN, SPDIAM, INFO )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       INTEGER            IFIRST, ILAST, INFO, N, OFFSET
012:       REAL               PIVMIN, RTOL, SPDIAM
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IWORK( * )
016:       REAL               D( * ), E2( * ), W( * ),
017:      $                   WERR( * ), WORK( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  Given the initial eigenvalue approximations of T, SLARRJ
024: *  does  bisection to refine the eigenvalues of T,
025: *  W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
026: *  guesses for these eigenvalues are input in W, the corresponding estimate
027: *  of the error in these guesses in WERR. During bisection, intervals
028: *  [left, right] are maintained by storing their mid-points and
029: *  semi-widths in the arrays W and WERR respectively.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  N       (input) INTEGER
035: *          The order of the matrix.
036: *
037: *  D       (input) REAL             array, dimension (N)
038: *          The N diagonal elements of T.
039: *
040: *  E2      (input) REAL             array, dimension (N-1)
041: *          The Squares of the (N-1) subdiagonal elements of T.
042: *
043: *  IFIRST  (input) INTEGER
044: *          The index of the first eigenvalue to be computed.
045: *
046: *  ILAST   (input) INTEGER
047: *          The index of the last eigenvalue to be computed.
048: *
049: *  RTOL   (input) REAL            
050: *          Tolerance for the convergence of the bisection intervals.
051: *          An interval [LEFT,RIGHT] has converged if
052: *          RIGHT-LEFT.LT.RTOL*MAX(|LEFT|,|RIGHT|).
053: *
054: *  OFFSET  (input) INTEGER
055: *          Offset for the arrays W and WERR, i.e., the IFIRST-OFFSET
056: *          through ILAST-OFFSET elements of these arrays are to be used.
057: *
058: *  W       (input/output) REAL             array, dimension (N)
059: *          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
060: *          estimates of the eigenvalues of L D L^T indexed IFIRST through
061: *          ILAST.
062: *          On output, these estimates are refined.
063: *
064: *  WERR    (input/output) REAL             array, dimension (N)
065: *          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
066: *          the errors in the estimates of the corresponding elements in W.
067: *          On output, these errors are refined.
068: *
069: *  WORK    (workspace) REAL             array, dimension (2*N)
070: *          Workspace.
071: *
072: *  IWORK   (workspace) INTEGER array, dimension (2*N)
073: *          Workspace.
074: *
075: *  PIVMIN  (input) DOUBLE PRECISION
076: *          The minimum pivot in the Sturm sequence for T.
077: *
078: *  SPDIAM  (input) DOUBLE PRECISION
079: *          The spectral diameter of T.
080: *
081: *  INFO    (output) INTEGER
082: *          Error flag.
083: *
084: *  Further Details
085: *  ===============
086: *
087: *  Based on contributions by
088: *     Beresford Parlett, University of California, Berkeley, USA
089: *     Jim Demmel, University of California, Berkeley, USA
090: *     Inderjit Dhillon, University of Texas, Austin, USA
091: *     Osni Marques, LBNL/NERSC, USA
092: *     Christof Voemel, University of California, Berkeley, USA
093: *
094: *  =====================================================================
095: *
096: *     .. Parameters ..
097:       REAL               ZERO, ONE, TWO, HALF
098:       PARAMETER        ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
099:      $                   HALF = 0.5E0 )
100:       INTEGER   MAXITR
101: *     ..
102: *     .. Local Scalars ..
103:       INTEGER            CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT,
104:      $                   OLNINT, P, PREV, SAVI1
105:       REAL               DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH
106: *
107: *     ..
108: *     .. Intrinsic Functions ..
109:       INTRINSIC          ABS, MAX
110: *     ..
111: *     .. Executable Statements ..
112: *
113:       INFO = 0
114: *
115:       MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
116:      $           LOG( TWO ) ) + 2
117: *
118: *     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
119: *     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
120: *     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
121: *     for an unconverged interval is set to the index of the next unconverged
122: *     interval, and is -1 or 0 for a converged interval. Thus a linked
123: *     list of unconverged intervals is set up.
124: *
125: 
126:       I1 = IFIRST
127:       I2 = ILAST
128: *     The number of unconverged intervals
129:       NINT = 0
130: *     The last unconverged interval found
131:       PREV = 0
132:       DO 75 I = I1, I2
133:          K = 2*I
134:          II = I - OFFSET
135:          LEFT = W( II ) - WERR( II )
136:          MID = W(II)
137:          RIGHT = W( II ) + WERR( II )
138:          WIDTH = RIGHT - MID
139:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
140: 
141: *        The following test prevents the test of converged intervals
142:          IF( WIDTH.LT.RTOL*TMP ) THEN
143: *           This interval has already converged and does not need refinement.
144: *           (Note that the gaps might change through refining the
145: *            eigenvalues, however, they can only get bigger.)
146: *           Remove it from the list.
147:             IWORK( K-1 ) = -1
148: *           Make sure that I1 always points to the first unconverged interval
149:             IF((I.EQ.I1).AND.(I.LT.I2)) I1 = I + 1
150:             IF((PREV.GE.I1).AND.(I.LE.I2)) IWORK( 2*PREV-1 ) = I + 1
151:          ELSE
152: *           unconverged interval found
153:             PREV = I
154: *           Make sure that [LEFT,RIGHT] contains the desired eigenvalue
155: *
156: *           Do while( CNT(LEFT).GT.I-1 )
157: *
158:             FAC = ONE
159:  20         CONTINUE
160:             CNT = 0
161:             S = LEFT
162:             DPLUS = D( 1 ) - S
163:             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
164:             DO 30 J = 2, N
165:                DPLUS = D( J ) - S - E2( J-1 )/DPLUS
166:                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
167:  30         CONTINUE
168:             IF( CNT.GT.I-1 ) THEN
169:                LEFT = LEFT - WERR( II )*FAC
170:                FAC = TWO*FAC
171:                GO TO 20
172:             END IF
173: *
174: *           Do while( CNT(RIGHT).LT.I )
175: *
176:             FAC = ONE
177:  50         CONTINUE
178:             CNT = 0
179:             S = RIGHT
180:             DPLUS = D( 1 ) - S
181:             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
182:             DO 60 J = 2, N
183:                DPLUS = D( J ) - S - E2( J-1 )/DPLUS
184:                IF( DPLUS.LT.ZERO ) CNT = CNT + 1
185:  60         CONTINUE
186:             IF( CNT.LT.I ) THEN
187:                RIGHT = RIGHT + WERR( II )*FAC
188:                FAC = TWO*FAC
189:                GO TO 50
190:             END IF
191:             NINT = NINT + 1
192:             IWORK( K-1 ) = I + 1
193:             IWORK( K ) = CNT
194:          END IF
195:          WORK( K-1 ) = LEFT
196:          WORK( K ) = RIGHT
197:  75   CONTINUE
198: 
199: 
200:       SAVI1 = I1
201: *
202: *     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
203: *     and while (ITER.LT.MAXITR)
204: *
205:       ITER = 0
206:  80   CONTINUE
207:       PREV = I1 - 1
208:       I = I1
209:       OLNINT = NINT
210: 
211:       DO 100 P = 1, OLNINT
212:          K = 2*I
213:          II = I - OFFSET
214:          NEXT = IWORK( K-1 )
215:          LEFT = WORK( K-1 )
216:          RIGHT = WORK( K )
217:          MID = HALF*( LEFT + RIGHT )
218: 
219: *        semiwidth of interval
220:          WIDTH = RIGHT - MID
221:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
222: 
223:          IF( ( WIDTH.LT.RTOL*TMP ) .OR.
224:      $      (ITER.EQ.MAXITR) )THEN
225: *           reduce number of unconverged intervals
226:             NINT = NINT - 1
227: *           Mark interval as converged.
228:             IWORK( K-1 ) = 0
229:             IF( I1.EQ.I ) THEN
230:                I1 = NEXT
231:             ELSE
232: *              Prev holds the last unconverged interval previously examined
233:                IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
234:             END IF
235:             I = NEXT
236:             GO TO 100
237:          END IF
238:          PREV = I
239: *
240: *        Perform one bisection step
241: *
242:          CNT = 0
243:          S = MID
244:          DPLUS = D( 1 ) - S
245:          IF( DPLUS.LT.ZERO ) CNT = CNT + 1
246:          DO 90 J = 2, N
247:             DPLUS = D( J ) - S - E2( J-1 )/DPLUS
248:             IF( DPLUS.LT.ZERO ) CNT = CNT + 1
249:  90      CONTINUE
250:          IF( CNT.LE.I-1 ) THEN
251:             WORK( K-1 ) = MID
252:          ELSE
253:             WORK( K ) = MID
254:          END IF
255:          I = NEXT
256: 
257:  100  CONTINUE
258:       ITER = ITER + 1
259: *     do another loop if there are still unconverged intervals
260: *     However, in the last iteration, all intervals are accepted
261: *     since this is the best we can do.
262:       IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
263: *
264: *
265: *     At this point, all the intervals have converged
266:       DO 110 I = SAVI1, ILAST
267:          K = 2*I
268:          II = I - OFFSET
269: *        All intervals marked by '0' have been refined.
270:          IF( IWORK( K-1 ).EQ.0 ) THEN
271:             W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
272:             WERR( II ) = WORK( K ) - W( II )
273:          END IF
274:  110  CONTINUE
275: *
276: 
277:       RETURN
278: *
279: *     End of SLARRJ
280: *
281:       END
282: