001:       SUBROUTINE ZLAR1V( 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:       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
013:      $                   RQCORR, ZTZ
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            ISUPPZ( * )
017:       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
018:      $                  WORK( * )
019:       COMPLEX*16       Z( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZLAR1V 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) DOUBLE PRECISION
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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
063: *           The n diagonal elements of the diagonal matrix D.
064: *
065: *  LD       (input) DOUBLE PRECISION array, dimension (N-1)
066: *           The n-1 elements L(i)*D(i).
067: *
068: *  LLD      (input) DOUBLE PRECISION array, dimension (N-1)
069: *           The n-1 elements L(i)*L(i)*D(i).
070: *
071: *  PIVMIN   (input) DOUBLE PRECISION
072: *           The minimum pivot in the Sturm sequence.
073: *
074: *  GAPTOL   (input) DOUBLE PRECISION
075: *           Tolerance that indicates when eigenvector entries are negligible
076: *           w.r.t. their contribution to the residual.
077: *
078: *  Z        (input/output) COMPLEX*16       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) DOUBLE PRECISION
091: *           The square of the 2-norm of Z.
092: *
093: *  MINGMA   (output) DOUBLE PRECISION
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) DOUBLE PRECISION
112: *           NRMINV = 1/SQRT( ZTZ )
113: *
114: *  RESID    (output) DOUBLE PRECISION
115: *           The residual of the FP vector.
116: *           RESID = ABS( MINGMA )/SQRT( ZTZ )
117: *
118: *  RQCORR   (output) DOUBLE PRECISION
119: *           The Rayleigh Quotient correction to LAMBDA.
120: *           RQCORR = MINGMA*TMP
121: *
122: *  WORK     (workspace) DOUBLE PRECISION 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:       DOUBLE PRECISION   ZERO, ONE
138:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
139:       COMPLEX*16         CONE
140:       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
141: 
142: *     ..
143: *     .. Local Scalars ..
144:       LOGICAL            SAWNAN1, SAWNAN2
145:       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
146:      $                   R2
147:       DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
148: *     ..
149: *     .. External Functions ..
150:       LOGICAL DISNAN
151:       DOUBLE PRECISION   DLAMCH
152:       EXTERNAL           DISNAN, DLAMCH
153: *     ..
154: *     .. Intrinsic Functions ..
155:       INTRINSIC          ABS, DBLE
156: *     ..
157: *     .. Executable Statements ..
158: *
159:       EPS = DLAMCH( 'Precision' )
160: 
161: 
162:       IF( R.EQ.0 ) THEN
163:          R1 = B1
164:          R2 = BN
165:       ELSE
166:          R1 = R
167:          R2 = R
168:       END IF
169: 
170: *     Storage for LPLUS
171:       INDLPL = 0
172: *     Storage for UMINUS
173:       INDUMN = N
174:       INDS = 2*N + 1
175:       INDP = 3*N + 1
176: 
177:       IF( B1.EQ.1 ) THEN
178:          WORK( INDS ) = ZERO
179:       ELSE
180:          WORK( INDS+B1-1 ) = LLD( B1-1 )
181:       END IF
182: 
183: *
184: *     Compute the stationary transform (using the differential form)
185: *     until the index R2.
186: *
187:       SAWNAN1 = .FALSE.
188:       NEG1 = 0
189:       S = WORK( INDS+B1-1 ) - LAMBDA
190:       DO 50 I = B1, R1 - 1
191:          DPLUS = D( I ) + S
192:          WORK( INDLPL+I ) = LD( I ) / DPLUS
193:          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
194:          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
195:          S = WORK( INDS+I ) - LAMBDA
196:  50   CONTINUE
197:       SAWNAN1 = DISNAN( S )
198:       IF( SAWNAN1 ) GOTO 60
199:       DO 51 I = R1, R2 - 1
200:          DPLUS = D( I ) + S
201:          WORK( INDLPL+I ) = LD( I ) / DPLUS
202:          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
203:          S = WORK( INDS+I ) - LAMBDA
204:  51   CONTINUE
205:       SAWNAN1 = DISNAN( S )
206: *
207:  60   CONTINUE
208:       IF( SAWNAN1 ) THEN
209: *        Runs a slower version of the above loop if a NaN is detected
210:          NEG1 = 0
211:          S = WORK( INDS+B1-1 ) - LAMBDA
212:          DO 70 I = B1, R1 - 1
213:             DPLUS = D( I ) + S
214:             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
215:             WORK( INDLPL+I ) = LD( I ) / DPLUS
216:             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
217:             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
218:             IF( WORK( INDLPL+I ).EQ.ZERO )
219:      $                      WORK( INDS+I ) = LLD( I )
220:             S = WORK( INDS+I ) - LAMBDA
221:  70      CONTINUE
222:          DO 71 I = R1, R2 - 1
223:             DPLUS = D( I ) + S
224:             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
225:             WORK( INDLPL+I ) = LD( I ) / DPLUS
226:             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
227:             IF( WORK( INDLPL+I ).EQ.ZERO )
228:      $                      WORK( INDS+I ) = LLD( I )
229:             S = WORK( INDS+I ) - LAMBDA
230:  71      CONTINUE
231:       END IF
232: *
233: *     Compute the progressive transform (using the differential form)
234: *     until the index R1
235: *
236:       SAWNAN2 = .FALSE.
237:       NEG2 = 0
238:       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
239:       DO 80 I = BN - 1, R1, -1
240:          DMINUS = LLD( I ) + WORK( INDP+I )
241:          TMP = D( I ) / DMINUS
242:          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
243:          WORK( INDUMN+I ) = L( I )*TMP
244:          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
245:  80   CONTINUE
246:       TMP = WORK( INDP+R1-1 )
247:       SAWNAN2 = DISNAN( TMP )
248: 
249:       IF( SAWNAN2 ) THEN
250: *        Runs a slower version of the above loop if a NaN is detected
251:          NEG2 = 0
252:          DO 100 I = BN-1, R1, -1
253:             DMINUS = LLD( I ) + WORK( INDP+I )
254:             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
255:             TMP = D( I ) / DMINUS
256:             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
257:             WORK( INDUMN+I ) = L( I )*TMP
258:             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
259:             IF( TMP.EQ.ZERO )
260:      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
261:  100     CONTINUE
262:       END IF
263: *
264: *     Find the index (from R1 to R2) of the largest (in magnitude)
265: *     diagonal element of the inverse
266: *
267:       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
268:       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
269:       IF( WANTNC ) THEN
270:          NEGCNT = NEG1 + NEG2
271:       ELSE
272:          NEGCNT = -1
273:       ENDIF
274:       IF( ABS(MINGMA).EQ.ZERO )
275:      $   MINGMA = EPS*WORK( INDS+R1-1 )
276:       R = R1
277:       DO 110 I = R1, R2 - 1
278:          TMP = WORK( INDS+I ) + WORK( INDP+I )
279:          IF( TMP.EQ.ZERO )
280:      $      TMP = EPS*WORK( INDS+I )
281:          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
282:             MINGMA = TMP
283:             R = I + 1
284:          END IF
285:  110  CONTINUE
286: *
287: *     Compute the FP vector: solve N^T v = e_r
288: *
289:       ISUPPZ( 1 ) = B1
290:       ISUPPZ( 2 ) = BN
291:       Z( R ) = CONE
292:       ZTZ = ONE
293: *
294: *     Compute the FP vector upwards from R
295: *
296:       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
297:          DO 210 I = R-1, B1, -1
298:             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
299:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
300:      $           THEN
301:                Z( I ) = ZERO
302:                ISUPPZ( 1 ) = I + 1
303:                GOTO 220
304:             ENDIF
305:             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
306:  210     CONTINUE
307:  220     CONTINUE
308:       ELSE
309: *        Run slower loop if NaN occurred.
310:          DO 230 I = R - 1, B1, -1
311:             IF( Z( I+1 ).EQ.ZERO ) THEN
312:                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
313:             ELSE
314:                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
315:             END IF
316:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
317:      $           THEN
318:                Z( I ) = ZERO
319:                ISUPPZ( 1 ) = I + 1
320:                GO TO 240
321:             END IF
322:             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
323:  230     CONTINUE
324:  240     CONTINUE
325:       ENDIF
326: 
327: *     Compute the FP vector downwards from R in blocks of size BLKSIZ
328:       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
329:          DO 250 I = R, BN-1
330:             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
331:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
332:      $         THEN
333:                Z( I+1 ) = ZERO
334:                ISUPPZ( 2 ) = I
335:                GO TO 260
336:             END IF
337:             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
338:  250     CONTINUE
339:  260     CONTINUE
340:       ELSE
341: *        Run slower loop if NaN occurred.
342:          DO 270 I = R, BN - 1
343:             IF( Z( I ).EQ.ZERO ) THEN
344:                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
345:             ELSE
346:                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
347:             END IF
348:             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
349:      $           THEN
350:                Z( I+1 ) = ZERO
351:                ISUPPZ( 2 ) = I
352:                GO TO 280
353:             END IF
354:             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
355:  270     CONTINUE
356:  280     CONTINUE
357:       END IF
358: *
359: *     Compute quantities for convergence test
360: *
361:       TMP = ONE / ZTZ
362:       NRMINV = SQRT( TMP )
363:       RESID = ABS( MINGMA )*NRMINV
364:       RQCORR = MINGMA*TMP
365: *
366: *
367:       RETURN
368: *
369: *     End of ZLAR1V
370: *
371:       END
372: