OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
slar1va.f
Go to the documentation of this file.
1 SUBROUTINE slar1va(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 IMPLICIT NONE
6*
7* -- ScaLAPACK computational routine (version 2.0) --
8* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
9* July 4, 2010
10*
11* .. Scalar Arguments ..
12 LOGICAL WANTNC
13 INTEGER B1, BN, N, NEGCNT, R
14 REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
15 $ rqcorr, ztz
16* ..
17* .. Array Arguments ..
18 INTEGER ISUPPZ( * )
19 REAL D( * ), L( * ), LD( * ), LLD( * ),
20 $ work( * )
21 REAL Z( * )
22*
23* Purpose
24* =======
25*
26* SLAR1VA 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) REAL
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) REAL 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) REAL array, dimension (N)
64* The n diagonal elements of the diagonal matrix D.
65*
66* LD (input) REAL array, dimension (N-1)
67* The n-1 elements L(i)*D(i).
68*
69* LLD (input) REAL array, dimension (N-1)
70* The n-1 elements L(i)*L(i)*D(i).
71*
72* PIVMIN (input) REAL
73* The minimum pivot in the Sturm sequence.
74*
75* GAPTOL (input) REAL
76* Tolerance that indicates when eigenvector entries are negligible
77* w.r.t. their contribution to the residual.
78*
79* Z (input/output) REAL 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) REAL
92* The square of the 2-norm of Z.
93*
94* MINGMA (output) REAL
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) 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 INTEGER BLKLEN
139 PARAMETER ( BLKLEN = 16 )
140 REAL ZERO, ONE
141 parameter( zero = 0.0e0, one = 1.0e0 )
142
143* ..
144* .. Local Scalars ..
145 LOGICAL SAWNAN1, SAWNAN2
146 INTEGER BI, I, INDLPL, INDP, INDS, INDUMN, NB, NEG1,
147 $ neg2, nx, r1, r2, to
148 REAL ABSZCUR, ABSZPREV, DMINUS, DPLUS, EPS,
149 $ S, TMP, ZPREV
150* ..
151* .. External Functions ..
152 LOGICAL SISNAN
153 REAL SLAMCH
154 EXTERNAL sisnan, slamch
155* ..
156* .. Intrinsic Functions ..
157 INTRINSIC abs, max, min, real
158* ..
159* .. Executable Statements ..
160*
161 eps = slamch( 'Precision' )
162
163
164 IF( r.EQ.0 ) THEN
165 r1 = b1
166 r2 = bn
167 ELSE
168 r1 = r
169 r2 = r
170 END IF
171
172* Storage for LPLUS
173 indlpl = 0
174* Storage for UMINUS
175 indumn = n
176 inds = 2*n + 1
177 indp = 3*n + 1
178
179 IF( b1.EQ.1 ) THEN
180 work( inds ) = zero
181 ELSE
182 work( inds+b1-1 ) = lld( b1-1 )
183 END IF
184
185*
186* Compute the stationary transform (using the differential form)
187* until the index R2.
188*
189 sawnan1 = .false.
190 neg1 = 0
191 s = work( inds+b1-1 ) - lambda
192 DO 50 i = b1, r1 - 1
193 dplus = d( i ) + s
194 work( indlpl+i ) = ld( i ) / dplus
195 IF(dplus.LT.zero) neg1 = neg1 + 1
196 work( inds+i ) = s*work( indlpl+i )*l( i )
197 s = work( inds+i ) - lambda
198 50 CONTINUE
199 sawnan1 = sisnan( s )
200 IF( sawnan1 ) GOTO 60
201 DO 51 i = r1, r2 - 1
202 dplus = d( i ) + s
203 work( indlpl+i ) = ld( i ) / dplus
204 work( inds+i ) = s*work( indlpl+i )*l( i )
205 s = work( inds+i ) - lambda
206 51 CONTINUE
207 sawnan1 = sisnan( s )
208*
209 60 CONTINUE
210 IF( sawnan1 ) THEN
211* Runs a slower version of the above loop if a NaN is detected
212 neg1 = 0
213 s = work( inds+b1-1 ) - lambda
214 DO 70 i = b1, r1 - 1
215 dplus = d( i ) + s
216 IF(abs(dplus).LT.pivmin) dplus = -pivmin
217 work( indlpl+i ) = ld( i ) / dplus
218 IF(dplus.LT.zero) neg1 = neg1 + 1
219 work( inds+i ) = s*work( indlpl+i )*l( i )
220 IF( work( indlpl+i ).EQ.zero )
221 $ work( inds+i ) = lld( i )
222 s = work( inds+i ) - lambda
223 70 CONTINUE
224 DO 71 i = r1, r2 - 1
225 dplus = d( i ) + s
226 IF(abs(dplus).LT.pivmin) dplus = -pivmin
227 work( indlpl+i ) = ld( i ) / dplus
228 work( inds+i ) = s*work( indlpl+i )*l( i )
229 IF( work( indlpl+i ).EQ.zero )
230 $ work( inds+i ) = lld( i )
231 s = work( inds+i ) - lambda
232 71 CONTINUE
233 END IF
234*
235* Compute the progressive transform (using the differential form)
236* until the index R1
237*
238 sawnan2 = .false.
239 neg2 = 0
240 work( indp+bn-1 ) = d( bn ) - lambda
241 DO 80 i = bn - 1, r1, -1
242 dminus = lld( i ) + work( indp+i )
243 tmp = d( i ) / dminus
244 IF(dminus.LT.zero) neg2 = neg2 + 1
245 work( indumn+i ) = l( i )*tmp
246 work( indp+i-1 ) = work( indp+i )*tmp - lambda
247 80 CONTINUE
248 tmp = work( indp+r1-1 )
249 sawnan2 = sisnan( tmp )
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 ) = one
293 ztz = one
294*
295* Compute the FP vector upwards from R
296*
297 nb = int((r-b1)/blklen)
298 nx = r-nb*blklen
299 IF( .NOT.sawnan1 ) THEN
300 DO 210 bi = r-1, nx, -blklen
301 to = bi-blklen+1
302 DO 205 i = bi, to, -1
303 z( i ) = -( work(indlpl+i)*z(i+1) )
304 ztz = ztz + z( i )*z( i )
305 205 CONTINUE
306 IF( abs(z(to)).LT.eps .AND.
307 $ abs(z(to+1)).LT.eps ) THEN
308 isuppz(1) = to
309 GOTO 220
310 ENDIF
311 210 CONTINUE
312 DO 215 i = nx-1, b1, -1
313 z( i ) = -( work(indlpl+i)*z(i+1) )
314 ztz = ztz + z( i )*z( i )
315 215 CONTINUE
316 220 CONTINUE
317 ELSE
318* Run slower loop if NaN occurred.
319 DO 230 bi = r-1, nx, -blklen
320 to = bi-blklen+1
321 DO 225 i = bi, to, -1
322 IF( z( i+1 ).EQ.zero ) THEN
323 z( i ) = -( ld( i+1 ) / ld( i ) )*z( i+2 )
324 ELSE
325 z( i ) = -( work( indlpl+i )*z( i+1 ) )
326 END IF
327 ztz = ztz + z( i )*z( i )
328 225 CONTINUE
329 IF( abs(z(to)).LT.eps .AND.
330 $ abs(z(to+1)).LT.eps ) THEN
331 isuppz(1) = to
332 GOTO 240
333 ENDIF
334 230 CONTINUE
335 DO 235 i = nx-1, b1, -1
336 IF( z( i+1 ).EQ.zero ) THEN
337 z( i ) = -( ld( i+1 ) / ld( i ) )*z( i+2 )
338 ELSE
339 z( i ) = -( work( indlpl+i )*z( i+1 ) )
340 END IF
341 ztz = ztz + z( i )*z( i )
342 235 CONTINUE
343 240 CONTINUE
344 ENDIF
345 DO 245 i= b1, (isuppz(1)-1)
346 z(i) = zero
347 245 CONTINUE
348
349* Compute the FP vector downwards from R in blocks of size BLKLEN
350 IF( .NOT.sawnan2 ) THEN
351 DO 260 bi = r+1, bn, blklen
352 to = bi+blklen-1
353 IF ( to.LE.bn ) THEN
354 DO 250 i = bi, to
355 z(i) = -(work(indumn+i-1)*z(i-1))
356 ztz = ztz + z( i )*z( i )
357 250 CONTINUE
358 IF( abs(z(to)).LE.eps .AND.
359 $ abs(z(to-1)).LE.eps ) THEN
360 isuppz(2) = to
361 GOTO 265
362 ENDIF
363 ELSE
364 DO 255 i = bi, bn
365 z(i) = -(work(indumn+i-1)*z(i-1))
366 ztz = ztz + z( i )*z( i )
367 255 CONTINUE
368 ENDIF
369 260 CONTINUE
370 265 CONTINUE
371 ELSE
372* Run slower loop if NaN occurred.
373 DO 280 bi = r+1, bn, blklen
374 to = bi+blklen-1
375 IF ( to.LE.bn ) THEN
376 DO 270 i = bi, to
377 zprev = z(i-1)
378 abszprev = abs(zprev)
379 IF( zprev.NE.zero ) THEN
380 z(i)= -(work(indumn+i-1)*zprev)
381 ELSE
382 z(i)= -(ld(i-2)/ld(i-1))*z(i-2)
383 END IF
384 abszcur = abs(z(i))
385 ztz = ztz + abszcur**2
386 270 CONTINUE
387 IF( abszcur.LT.eps .AND.
388 $ abszprev.LT.eps ) THEN
389 isuppz(2) = i
390 GOTO 285
391 ENDIF
392 ELSE
393 DO 275 i = bi, bn
394 zprev = z(i-1)
395 abszprev = abs(zprev)
396 IF( zprev.NE.zero ) THEN
397 z(i)= -(work(indumn+i-1)*zprev)
398 ELSE
399 z(i)= -(ld(i-2)/ld(i-1))*z(i-2)
400 END IF
401 abszcur = abs(z(i))
402 ztz = ztz + abszcur**2
403 275 CONTINUE
404 ENDIF
405 280 CONTINUE
406 285 CONTINUE
407 END IF
408 DO 290 i= isuppz(2)+1,bn
409 z(i) = zero
410 290 CONTINUE
411*
412* Compute quantities for convergence test
413*
414 tmp = one / ztz
415 nrminv = sqrt( tmp )
416 resid = abs( mingma )*nrminv
417 rqcorr = mingma*tmp
418*
419 RETURN
420*
421* End of SLAR1VA
422*
423 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
slar1va
subroutine slar1va(n, b1, bn, lambda, d, l, ld, lld, pivmin, gaptol, z, wantnc, negcnt, ztz, mingma, r, isuppz, nrminv, resid, rqcorr, work)
Definition slar1va.f:4