1 SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
2 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
3 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 LOGICAL WANTNC
12 INTEGER B1, BN, N, NEGCNT, R
13 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
14 $ RQCORR, ZTZ
15 * ..
16 * .. Array Arguments ..
17 INTEGER ISUPPZ( * )
18 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
19 $ WORK( * )
20 COMPLEX*16 Z( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZLAR1V computes the (scaled) r-th column of the inverse of
27 * the sumbmatrix in rows B1 through BN of the tridiagonal matrix
28 * L D L**T - sigma I. When sigma is close to an eigenvalue, the
29 * computed vector is an accurate eigenvector. Usually, r corresponds
30 * to the index where the eigenvector is largest in magnitude.
31 * The following steps accomplish this computation :
32 * (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
33 * (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
34 * (c) Computation of the diagonal elements of the inverse of
35 * L D L**T - sigma I by combining the above transforms, and choosing
36 * r as the index where the diagonal of the inverse is (one of the)
37 * largest in magnitude.
38 * (d) Computation of the (scaled) r-th column of the inverse using the
39 * twisted factorization obtained by combining the top part of the
40 * the stationary and the bottom part of the progressive transform.
41 *
42 * Arguments
43 * =========
44 *
45 * N (input) INTEGER
46 * The order of the matrix L D L**T.
47 *
48 * B1 (input) INTEGER
49 * First index of the submatrix of L D L**T.
50 *
51 * BN (input) INTEGER
52 * Last index of the submatrix of L D L**T.
53 *
54 * LAMBDA (input) DOUBLE PRECISION
55 * The shift. In order to compute an accurate eigenvector,
56 * LAMBDA should be a good approximation to an eigenvalue
57 * of L D L**T.
58 *
59 * L (input) DOUBLE PRECISION array, dimension (N-1)
60 * The (n-1) subdiagonal elements of the unit bidiagonal matrix
61 * L, in elements 1 to N-1.
62 *
63 * D (input) DOUBLE PRECISION array, dimension (N)
64 * The n diagonal elements of the diagonal matrix D.
65 *
66 * LD (input) DOUBLE PRECISION array, dimension (N-1)
67 * The n-1 elements L(i)*D(i).
68 *
69 * LLD (input) DOUBLE PRECISION array, dimension (N-1)
70 * The n-1 elements L(i)*L(i)*D(i).
71 *
72 * PIVMIN (input) DOUBLE PRECISION
73 * The minimum pivot in the Sturm sequence.
74 *
75 * GAPTOL (input) DOUBLE PRECISION
76 * Tolerance that indicates when eigenvector entries are negligible
77 * w.r.t. their contribution to the residual.
78 *
79 * Z (input/output) COMPLEX*16 array, dimension (N)
80 * On input, all entries of Z must be set to 0.
81 * On output, Z contains the (scaled) r-th column of the
82 * inverse. The scaling is such that Z(R) equals 1.
83 *
84 * WANTNC (input) LOGICAL
85 * Specifies whether NEGCNT has to be computed.
86 *
87 * NEGCNT (output) INTEGER
88 * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
89 * in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
90 *
91 * ZTZ (output) DOUBLE PRECISION
92 * The square of the 2-norm of Z.
93 *
94 * MINGMA (output) DOUBLE PRECISION
95 * The reciprocal of the largest (in magnitude) diagonal
96 * element of the inverse of L D L**T - sigma I.
97 *
98 * R (input/output) INTEGER
99 * 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
2 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
3 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 LOGICAL WANTNC
12 INTEGER B1, BN, N, NEGCNT, R
13 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
14 $ RQCORR, ZTZ
15 * ..
16 * .. Array Arguments ..
17 INTEGER ISUPPZ( * )
18 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
19 $ WORK( * )
20 COMPLEX*16 Z( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * ZLAR1V computes the (scaled) r-th column of the inverse of
27 * the sumbmatrix in rows B1 through BN of the tridiagonal matrix
28 * L D L**T - sigma I. When sigma is close to an eigenvalue, the
29 * computed vector is an accurate eigenvector. Usually, r corresponds
30 * to the index where the eigenvector is largest in magnitude.
31 * The following steps accomplish this computation :
32 * (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
33 * (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
34 * (c) Computation of the diagonal elements of the inverse of
35 * L D L**T - sigma I by combining the above transforms, and choosing
36 * r as the index where the diagonal of the inverse is (one of the)
37 * largest in magnitude.
38 * (d) Computation of the (scaled) r-th column of the inverse using the
39 * twisted factorization obtained by combining the top part of the
40 * the stationary and the bottom part of the progressive transform.
41 *
42 * Arguments
43 * =========
44 *
45 * N (input) INTEGER
46 * The order of the matrix L D L**T.
47 *
48 * B1 (input) INTEGER
49 * First index of the submatrix of L D L**T.
50 *
51 * BN (input) INTEGER
52 * Last index of the submatrix of L D L**T.
53 *
54 * LAMBDA (input) DOUBLE PRECISION
55 * The shift. In order to compute an accurate eigenvector,
56 * LAMBDA should be a good approximation to an eigenvalue
57 * of L D L**T.
58 *
59 * L (input) DOUBLE PRECISION array, dimension (N-1)
60 * The (n-1) subdiagonal elements of the unit bidiagonal matrix
61 * L, in elements 1 to N-1.
62 *
63 * D (input) DOUBLE PRECISION array, dimension (N)
64 * The n diagonal elements of the diagonal matrix D.
65 *
66 * LD (input) DOUBLE PRECISION array, dimension (N-1)
67 * The n-1 elements L(i)*D(i).
68 *
69 * LLD (input) DOUBLE PRECISION array, dimension (N-1)
70 * The n-1 elements L(i)*L(i)*D(i).
71 *
72 * PIVMIN (input) DOUBLE PRECISION
73 * The minimum pivot in the Sturm sequence.
74 *
75 * GAPTOL (input) DOUBLE PRECISION
76 * Tolerance that indicates when eigenvector entries are negligible
77 * w.r.t. their contribution to the residual.
78 *
79 * Z (input/output) COMPLEX*16 array, dimension (N)
80 * On input, all entries of Z must be set to 0.
81 * On output, Z contains the (scaled) r-th column of the
82 * inverse. The scaling is such that Z(R) equals 1.
83 *
84 * WANTNC (input) LOGICAL
85 * Specifies whether NEGCNT has to be computed.
86 *
87 * NEGCNT (output) INTEGER
88 * If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
89 * in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
90 *
91 * ZTZ (output) DOUBLE PRECISION
92 * The square of the 2-norm of Z.
93 *
94 * MINGMA (output) DOUBLE PRECISION
95 * The reciprocal of the largest (in magnitude) diagonal
96 * element of the inverse of L D L**T - sigma I.
97 *
98 * R (input/output) INTEGER
99 * 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