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