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