001:       SUBROUTINE SLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
002:      $           PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
003:      $           R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            WANTNC
011:       INTEGER   B1, BN, N, NEGCNT, R
012:       REAL               GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
013:      $                   RQCORR, ZTZ
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            ISUPPZ( * )
017:       REAL               D( * ), L( * ), LD( * ), LLD( * ),
018:      $                  WORK( * )
019:       REAL             Z( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  SLAR1V computes the (scaled) r-th column of the inverse of
026: *  the sumbmatrix in rows B1 through BN of the tridiagonal matrix
027: *  L D L^T - sigma I. When sigma is close to an eigenvalue, the
028: *  computed vector is an accurate eigenvector. Usually, r corresponds
029: *  to the index where the eigenvector is largest in magnitude.
030: *  The following steps accomplish this computation :
031: *  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T,
032: *  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
033: *  (c) Computation of the diagonal elements of the inverse of
034: *      L D L^T - sigma I by combining the above transforms, and choosing
035: *      r as the index where the diagonal of the inverse is (one of the)
036: *      largest in magnitude.
037: *  (d) Computation of the (scaled) r-th column of the inverse using the
038: *      twisted factorization obtained by combining the top part of the
039: *      the stationary and the bottom part of the progressive transform.
040: *
041: *  Arguments
042: *  =========
043: *
044: *  N        (input) INTEGER
045: *           The order of the matrix L D L^T.
046: *
047: *  B1       (input) INTEGER
048: *           First index of the submatrix of L D L^T.
049: *
050: *  BN       (input) INTEGER
051: *           Last index of the submatrix of L D L^T.
052: *
053: *  LAMBDA    (input) REAL            
054: *           The shift. In order to compute an accurate eigenvector,
055: *           LAMBDA should be a good approximation to an eigenvalue
056: *           of L D L^T.
057: *
058: *  L        (input) REAL             array, dimension (N-1)
059: *           The (n-1) subdiagonal elements of the unit bidiagonal matrix
060: *           L, in elements 1 to N-1.
061: *
062: *  D        (input) REAL             array, dimension (N)
063: *           The n diagonal elements of the diagonal matrix D.
064: *
065: *  LD       (input) REAL             array, dimension (N-1)
066: *           The n-1 elements L(i)*D(i).
067: *
068: *  LLD      (input) REAL             array, dimension (N-1)
069: *           The n-1 elements L(i)*L(i)*D(i).
070: *
071: *  PIVMIN   (input) REAL            
072: *           The minimum pivot in the Sturm sequence.
073: *
074: *  GAPTOL   (input) REAL            
075: *           Tolerance that indicates when eigenvector entries are negligible
076: *           w.r.t. their contribution to the residual.
077: *
078: *  Z        (input/output) REAL             array, dimension (N)
079: *           On input, all entries of Z must be set to 0.
080: *           On output, Z contains the (scaled) r-th column of the
081: *           inverse. The scaling is such that Z(R) equals 1.
082: *
083: *  WANTNC   (input) LOGICAL
084: *           Specifies whether NEGCNT has to be computed.
085: *
086: *  NEGCNT   (output) INTEGER
087: *           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
088: *           in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise.
089: *
090: *  ZTZ      (output) REAL            
091: *           The square of the 2-norm of Z.
092: *
093: *  MINGMA   (output) REAL            
094: *           The reciprocal of the largest (in magnitude) diagonal
095: *           element of the inverse of L D L^T - sigma I.
096: *
097: *  R        (input/output) INTEGER
098: *           The twist index for the twisted factorization used to
099: *           compute Z.
100: *           On input, 0 <= R <= N. If R is input as 0, R is set to
101: *           the index where (L D L^T - sigma I)^{-1} is largest
102: *           in magnitude. If 1 <= R <= N, R is unchanged.
103: *           On output, R contains the twist index used to compute Z.
104: *           Ideally, R designates the position of the maximum entry in the
105: *           eigenvector.
106: *
107: *  ISUPPZ   (output) INTEGER array, dimension (2)
108: *           The support of the vector in Z, i.e., the vector Z is
109: *           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
110: *
111: *  NRMINV   (output) REAL            
112: *           NRMINV = 1/SQRT( ZTZ )
113: *
114: *  RESID    (output) REAL            
115: *           The residual of the FP vector.
116: *           RESID = ABS( MINGMA )/SQRT( ZTZ )
117: *
118: *  RQCORR   (output) REAL            
119: *           The Rayleigh Quotient correction to LAMBDA.
120: *           RQCORR = MINGMA*TMP
121: *
122: *  WORK     (workspace) REAL             array, dimension (4*N)
123: *
124: *  Further Details
125: *  ===============
126: *
127: *  Based on contributions by
128: *     Beresford Parlett, University of California, Berkeley, USA
129: *     Jim Demmel, University of California, Berkeley, USA
130: *     Inderjit Dhillon, University of Texas, Austin, USA
131: *     Osni Marques, LBNL/NERSC, USA
132: *     Christof Voemel, University of California, Berkeley, USA
133: *
134: *  =====================================================================
135: *
136: *     .. Parameters ..
137:       REAL               ZERO, ONE
138:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
139: 
140: *     ..
141: *     .. Local Scalars ..
142:       LOGICAL            SAWNAN1, SAWNAN2
143:       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
144:      $                   R2
145:       REAL               DMINUS, DPLUS, EPS, S, TMP
146: *     ..
147: *     .. External Functions ..
148:       LOGICAL SISNAN
149:       REAL               SLAMCH
150:       EXTERNAL           SISNAN, SLAMCH
151: *     ..
152: *     .. Intrinsic Functions ..
153:       INTRINSIC          ABS
154: *     ..
155: *     .. Executable Statements ..
156: *
157:       EPS = SLAMCH( 'Precision' )
158: 
159: 
160:       IF( R.EQ.0 ) THEN
161:          R1 = B1
162:          R2 = BN
163:       ELSE
164:          R1 = R
165:          R2 = R
166:       END IF
167: 
168: *     Storage for LPLUS
169:       INDLPL = 0
170: *     Storage for UMINUS
171:       INDUMN = N
172:       INDS = 2*N + 1
173:       INDP = 3*N + 1
174: 
175:       IF( B1.EQ.1 ) THEN
176:          WORK( INDS ) = ZERO
177:       ELSE
178:          WORK( INDS+B1-1 ) = LLD( B1-1 )
179:       END IF
180: 
181: *
182: *     Compute the stationary transform (using the differential form)
183: *     until the index R2.
184: *
185:       SAWNAN1 = .FALSE.
186:       NEG1 = 0
187:       S = WORK( INDS+B1-1 ) - LAMBDA
188:       DO 50 I = B1, R1 - 1
189:          DPLUS = D( I ) + S
190:          WORK( INDLPL+I ) = LD( I ) / DPLUS
191:          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
192:          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
193:          S = WORK( INDS+I ) - LAMBDA
194:  50   CONTINUE
195:       SAWNAN1 = SISNAN( S )
196:       IF( SAWNAN1 ) GOTO 60
197:       DO 51 I = R1, R2 - 1
198:          DPLUS = D( I ) + S
199:          WORK( INDLPL+I ) = LD( I ) / DPLUS
200:          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
201:          S = WORK( INDS+I ) - LAMBDA
202:  51   CONTINUE
203:       SAWNAN1 = SISNAN( S )
204: *
205:  60   CONTINUE
206:       IF( SAWNAN1 ) THEN
207: *        Runs a slower version of the above loop if a NaN is detected
208:          NEG1 = 0
209:          S = WORK( INDS+B1-1 ) - LAMBDA
210:          DO 70 I = B1, R1 - 1
211:             DPLUS = D( I ) + S
212:             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
213:             WORK( INDLPL+I ) = LD( I ) / DPLUS
214:             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
215:             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
216:             IF( WORK( INDLPL+I ).EQ.ZERO )
217:      $                      WORK( INDS+I ) = LLD( I )
218:             S = WORK( INDS+I ) - LAMBDA
219:  70      CONTINUE
220:          DO 71 I = R1, R2 - 1
221:             DPLUS = D( I ) + S
222:             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
223:             WORK( INDLPL+I ) = LD( I ) / DPLUS
224:             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
225:             IF( WORK( INDLPL+I ).EQ.ZERO )
226:      $                      WORK( INDS+I ) = LLD( I )
227:             S = WORK( INDS+I ) - LAMBDA
228:  71      CONTINUE
229:       END IF
230: *
231: *     Compute the progressive transform (using the differential form)
232: *     until the index R1
233: *
234:       SAWNAN2 = .FALSE.
235:       NEG2 = 0
236:       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
237:       DO 80 I = BN - 1, R1, -1
238:          DMINUS = LLD( I ) + WORK( INDP+I )
239:          TMP = D( I ) / DMINUS
240:          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
241:          WORK( INDUMN+I ) = L( I )*TMP
242:          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
243:  80   CONTINUE
244:       TMP = WORK( INDP+R1-1 )
245:       SAWNAN2 = SISNAN( TMP )
246: 
247:       IF( SAWNAN2 ) THEN
248: *        Runs a slower version of the above loop if a NaN is detected
249:          NEG2 = 0
250:          DO 100 I = BN-1, R1, -1
251:             DMINUS = LLD( I ) + WORK( INDP+I )
252:             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
253:             TMP = D( I ) / DMINUS
254:             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
255:             WORK( INDUMN+I ) = L( I )*TMP
256:             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
257:             IF( TMP.EQ.ZERO )
258:      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
259:  100     CONTINUE
260:       END IF
261: *
262: *     Find the index (from R1 to R2) of the largest (in magnitude)
263: *     diagonal element of the inverse
264: *
265:       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
266:       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
267:       IF( WANTNC ) THEN
268:          NEGCNT = NEG1 + NEG2
269:       ELSE
270:          NEGCNT = -1
271:       ENDIF
272:       IF( ABS(MINGMA).EQ.ZERO )
273:      $   MINGMA = EPS*WORK( INDS+R1-1 )
274:       R = R1
275:       DO 110 I = R1, R2 - 1
276:          TMP = WORK( INDS+I ) + WORK( INDP+I )
277:          IF( TMP.EQ.ZERO )
278:      $      TMP = EPS*WORK( INDS+I )
279:          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
280:             MINGMA = TMP
281:             R = I + 1
282:          END IF
283:  110  CONTINUE
284: *
285: *     Compute the FP vector: solve N^T v = e_r
286: *
287:       ISUPPZ( 1 ) = B1
288:       ISUPPZ( 2 ) = BN
289:       Z( R ) = ONE
290:       ZTZ = ONE
291: *
292: *     Compute the FP vector upwards from R
293: *
294:       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
295:          DO 210 I = R-1, B1, -1
296:             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
297:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
298:      $           THEN
299:                Z( I ) = ZERO
300:                ISUPPZ( 1 ) = I + 1
301:                GOTO 220
302:             ENDIF
303:             ZTZ = ZTZ + Z( I )*Z( I )
304:  210     CONTINUE
305:  220     CONTINUE
306:       ELSE
307: *        Run slower loop if NaN occurred.
308:          DO 230 I = R - 1, B1, -1
309:             IF( Z( I+1 ).EQ.ZERO ) THEN
310:                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
311:             ELSE
312:                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
313:             END IF
314:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
315:      $           THEN
316:                Z( I ) = ZERO
317:                ISUPPZ( 1 ) = I + 1
318:                GO TO 240
319:             END IF
320:             ZTZ = ZTZ + Z( I )*Z( I )
321:  230     CONTINUE
322:  240     CONTINUE
323:       ENDIF
324: 
325: *     Compute the FP vector downwards from R in blocks of size BLKSIZ
326:       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
327:          DO 250 I = R, BN-1
328:             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
329:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
330:      $         THEN
331:                Z( I+1 ) = ZERO
332:                ISUPPZ( 2 ) = I
333:                GO TO 260
334:             END IF
335:             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
336:  250     CONTINUE
337:  260     CONTINUE
338:       ELSE
339: *        Run slower loop if NaN occurred.
340:          DO 270 I = R, BN - 1
341:             IF( Z( I ).EQ.ZERO ) THEN
342:                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
343:             ELSE
344:                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
345:             END IF
346:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
347:      $           THEN
348:                Z( I+1 ) = ZERO
349:                ISUPPZ( 2 ) = I
350:                GO TO 280
351:             END IF
352:             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
353:  270     CONTINUE
354:  280     CONTINUE
355:       END IF
356: *
357: *     Compute quantities for convergence test
358: *
359:       TMP = ONE / ZTZ
360:       NRMINV = SQRT( TMP )
361:       RESID = ABS( MINGMA )*NRMINV
362:       RQCORR = MINGMA*TMP
363: *
364: *
365:       RETURN
366: *
367: *     End of SLAR1V
368: *
369:       END
370: