LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clar1v.f
Go to the documentation of this file.
1*> \brief \b CLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CLAR1V + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clar1v.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clar1v.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clar1v.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
22* PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
23* R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
24*
25* .. Scalar Arguments ..
26* LOGICAL WANTNC
27* INTEGER B1, BN, N, NEGCNT, R
28* REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
29* $ RQCORR, ZTZ
30* ..
31* .. Array Arguments ..
32* INTEGER ISUPPZ( * )
33* REAL D( * ), L( * ), LD( * ), LLD( * ),
34* $ WORK( * )
35* COMPLEX Z( * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> CLAR1V computes the (scaled) r-th column of the inverse of
45*> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
46*> L D L**T - sigma I. When sigma is close to an eigenvalue, the
47*> computed vector is an accurate eigenvector. Usually, r corresponds
48*> to the index where the eigenvector is largest in magnitude.
49*> The following steps accomplish this computation :
50*> (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
51*> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
52*> (c) Computation of the diagonal elements of the inverse of
53*> L D L**T - sigma I by combining the above transforms, and choosing
54*> r as the index where the diagonal of the inverse is (one of the)
55*> largest in magnitude.
56*> (d) Computation of the (scaled) r-th column of the inverse using the
57*> twisted factorization obtained by combining the top part of the
58*> the stationary and the bottom part of the progressive transform.
59*> \endverbatim
60*
61* Arguments:
62* ==========
63*
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The order of the matrix L D L**T.
68*> \endverbatim
69*>
70*> \param[in] B1
71*> \verbatim
72*> B1 is INTEGER
73*> First index of the submatrix of L D L**T.
74*> \endverbatim
75*>
76*> \param[in] BN
77*> \verbatim
78*> BN is INTEGER
79*> Last index of the submatrix of L D L**T.
80*> \endverbatim
81*>
82*> \param[in] LAMBDA
83*> \verbatim
84*> LAMBDA is REAL
85*> The shift. In order to compute an accurate eigenvector,
86*> LAMBDA should be a good approximation to an eigenvalue
87*> of L D L**T.
88*> \endverbatim
89*>
90*> \param[in] L
91*> \verbatim
92*> L is REAL array, dimension (N-1)
93*> The (n-1) subdiagonal elements of the unit bidiagonal matrix
94*> L, in elements 1 to N-1.
95*> \endverbatim
96*>
97*> \param[in] D
98*> \verbatim
99*> D is REAL array, dimension (N)
100*> The n diagonal elements of the diagonal matrix D.
101*> \endverbatim
102*>
103*> \param[in] LD
104*> \verbatim
105*> LD is REAL array, dimension (N-1)
106*> The n-1 elements L(i)*D(i).
107*> \endverbatim
108*>
109*> \param[in] LLD
110*> \verbatim
111*> LLD is REAL array, dimension (N-1)
112*> The n-1 elements L(i)*L(i)*D(i).
113*> \endverbatim
114*>
115*> \param[in] PIVMIN
116*> \verbatim
117*> PIVMIN is REAL
118*> The minimum pivot in the Sturm sequence.
119*> \endverbatim
120*>
121*> \param[in] GAPTOL
122*> \verbatim
123*> GAPTOL is REAL
124*> Tolerance that indicates when eigenvector entries are negligible
125*> w.r.t. their contribution to the residual.
126*> \endverbatim
127*>
128*> \param[in,out] Z
129*> \verbatim
130*> Z is COMPLEX array, dimension (N)
131*> On input, all entries of Z must be set to 0.
132*> On output, Z contains the (scaled) r-th column of the
133*> inverse. The scaling is such that Z(R) equals 1.
134*> \endverbatim
135*>
136*> \param[in] WANTNC
137*> \verbatim
138*> WANTNC is LOGICAL
139*> Specifies whether NEGCNT has to be computed.
140*> \endverbatim
141*>
142*> \param[out] NEGCNT
143*> \verbatim
144*> NEGCNT is INTEGER
145*> If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
146*> in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
147*> \endverbatim
148*>
149*> \param[out] ZTZ
150*> \verbatim
151*> ZTZ is REAL
152*> The square of the 2-norm of Z.
153*> \endverbatim
154*>
155*> \param[out] MINGMA
156*> \verbatim
157*> MINGMA is REAL
158*> The reciprocal of the largest (in magnitude) diagonal
159*> element of the inverse of L D L**T - sigma I.
160*> \endverbatim
161*>
162*> \param[in,out] R
163*> \verbatim
164*> R is INTEGER
165*> The twist index for the twisted factorization used to
166*> compute Z.
167*> On input, 0 <= R <= N. If R is input as 0, R is set to
168*> the index where (L D L**T - sigma I)^{-1} is largest
169*> in magnitude. If 1 <= R <= N, R is unchanged.
170*> On output, R contains the twist index used to compute Z.
171*> Ideally, R designates the position of the maximum entry in the
172*> eigenvector.
173*> \endverbatim
174*>
175*> \param[out] ISUPPZ
176*> \verbatim
177*> ISUPPZ is INTEGER array, dimension (2)
178*> The support of the vector in Z, i.e., the vector Z is
179*> nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
180*> \endverbatim
181*>
182*> \param[out] NRMINV
183*> \verbatim
184*> NRMINV is REAL
185*> NRMINV = 1/SQRT( ZTZ )
186*> \endverbatim
187*>
188*> \param[out] RESID
189*> \verbatim
190*> RESID is REAL
191*> The residual of the FP vector.
192*> RESID = ABS( MINGMA )/SQRT( ZTZ )
193*> \endverbatim
194*>
195*> \param[out] RQCORR
196*> \verbatim
197*> RQCORR is REAL
198*> The Rayleigh Quotient correction to LAMBDA.
199*> RQCORR = MINGMA*TMP
200*> \endverbatim
201*>
202*> \param[out] WORK
203*> \verbatim
204*> WORK is REAL array, dimension (4*N)
205*> \endverbatim
206*
207* Authors:
208* ========
209*
210*> \author Univ. of Tennessee
211*> \author Univ. of California Berkeley
212*> \author Univ. of Colorado Denver
213*> \author NAG Ltd.
214*
215*> \ingroup lar1v
216*
217*> \par Contributors:
218* ==================
219*>
220*> Beresford Parlett, University of California, Berkeley, USA \n
221*> Jim Demmel, University of California, Berkeley, USA \n
222*> Inderjit Dhillon, University of Texas, Austin, USA \n
223*> Osni Marques, LBNL/NERSC, USA \n
224*> Christof Voemel, University of California, Berkeley, USA
225*
226* =====================================================================
227 SUBROUTINE clar1v( N, B1, BN, LAMBDA, D, L, LD, LLD,
228 $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
229 $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
230*
231* -- LAPACK auxiliary routine --
232* -- LAPACK is a software package provided by Univ. of Tennessee, --
233* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234*
235* .. Scalar Arguments ..
236 LOGICAL WANTNC
237 INTEGER B1, BN, N, NEGCNT, R
238 REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
239 $ rqcorr, ztz
240* ..
241* .. Array Arguments ..
242 INTEGER ISUPPZ( * )
243 REAL D( * ), L( * ), LD( * ), LLD( * ),
244 $ work( * )
245 COMPLEX Z( * )
246* ..
247*
248* =====================================================================
249*
250* .. Parameters ..
251 REAL ZERO, ONE
252 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
253 COMPLEX CONE
254 parameter( cone = ( 1.0e0, 0.0e0 ) )
255
256* ..
257* .. Local Scalars ..
258 LOGICAL SAWNAN1, SAWNAN2
259 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
260 $ r2
261 REAL DMINUS, DPLUS, EPS, S, TMP
262* ..
263* .. External Functions ..
264 LOGICAL SISNAN
265 REAL SLAMCH
266 EXTERNAL sisnan, slamch
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC abs, real
270* ..
271* .. Executable Statements ..
272*
273 eps = slamch( 'Precision' )
274
275
276 IF( r.EQ.0 ) THEN
277 r1 = b1
278 r2 = bn
279 ELSE
280 r1 = r
281 r2 = r
282 END IF
283
284* Storage for LPLUS
285 indlpl = 0
286* Storage for UMINUS
287 indumn = n
288 inds = 2*n + 1
289 indp = 3*n + 1
290
291 IF( b1.EQ.1 ) THEN
292 work( inds ) = zero
293 ELSE
294 work( inds+b1-1 ) = lld( b1-1 )
295 END IF
296
297*
298* Compute the stationary transform (using the differential form)
299* until the index R2.
300*
301 sawnan1 = .false.
302 neg1 = 0
303 s = work( inds+b1-1 ) - lambda
304 DO 50 i = b1, r1 - 1
305 dplus = d( i ) + s
306 work( indlpl+i ) = ld( i ) / dplus
307 IF(dplus.LT.zero) neg1 = neg1 + 1
308 work( inds+i ) = s*work( indlpl+i )*l( i )
309 s = work( inds+i ) - lambda
310 50 CONTINUE
311 sawnan1 = sisnan( s )
312 IF( sawnan1 ) GOTO 60
313 DO 51 i = r1, r2 - 1
314 dplus = d( i ) + s
315 work( indlpl+i ) = ld( i ) / dplus
316 work( inds+i ) = s*work( indlpl+i )*l( i )
317 s = work( inds+i ) - lambda
318 51 CONTINUE
319 sawnan1 = sisnan( s )
320*
321 60 CONTINUE
322 IF( sawnan1 ) THEN
323* Runs a slower version of the above loop if a NaN is detected
324 neg1 = 0
325 s = work( inds+b1-1 ) - lambda
326 DO 70 i = b1, r1 - 1
327 dplus = d( i ) + s
328 IF(abs(dplus).LT.pivmin) dplus = -pivmin
329 work( indlpl+i ) = ld( i ) / dplus
330 IF(dplus.LT.zero) neg1 = neg1 + 1
331 work( inds+i ) = s*work( indlpl+i )*l( i )
332 IF( work( indlpl+i ).EQ.zero )
333 $ work( inds+i ) = lld( i )
334 s = work( inds+i ) - lambda
335 70 CONTINUE
336 DO 71 i = r1, r2 - 1
337 dplus = d( i ) + s
338 IF(abs(dplus).LT.pivmin) dplus = -pivmin
339 work( indlpl+i ) = ld( i ) / dplus
340 work( inds+i ) = s*work( indlpl+i )*l( i )
341 IF( work( indlpl+i ).EQ.zero )
342 $ work( inds+i ) = lld( i )
343 s = work( inds+i ) - lambda
344 71 CONTINUE
345 END IF
346*
347* Compute the progressive transform (using the differential form)
348* until the index R1
349*
350 sawnan2 = .false.
351 neg2 = 0
352 work( indp+bn-1 ) = d( bn ) - lambda
353 DO 80 i = bn - 1, r1, -1
354 dminus = lld( i ) + work( indp+i )
355 tmp = d( i ) / dminus
356 IF(dminus.LT.zero) neg2 = neg2 + 1
357 work( indumn+i ) = l( i )*tmp
358 work( indp+i-1 ) = work( indp+i )*tmp - lambda
359 80 CONTINUE
360 tmp = work( indp+r1-1 )
361 sawnan2 = sisnan( tmp )
362
363 IF( sawnan2 ) THEN
364* Runs a slower version of the above loop if a NaN is detected
365 neg2 = 0
366 DO 100 i = bn-1, r1, -1
367 dminus = lld( i ) + work( indp+i )
368 IF(abs(dminus).LT.pivmin) dminus = -pivmin
369 tmp = d( i ) / dminus
370 IF(dminus.LT.zero) neg2 = neg2 + 1
371 work( indumn+i ) = l( i )*tmp
372 work( indp+i-1 ) = work( indp+i )*tmp - lambda
373 IF( tmp.EQ.zero )
374 $ work( indp+i-1 ) = d( i ) - lambda
375 100 CONTINUE
376 END IF
377*
378* Find the index (from R1 to R2) of the largest (in magnitude)
379* diagonal element of the inverse
380*
381 mingma = work( inds+r1-1 ) + work( indp+r1-1 )
382 IF( mingma.LT.zero ) neg1 = neg1 + 1
383 IF( wantnc ) THEN
384 negcnt = neg1 + neg2
385 ELSE
386 negcnt = -1
387 ENDIF
388 IF( abs(mingma).EQ.zero )
389 $ mingma = eps*work( inds+r1-1 )
390 r = r1
391 DO 110 i = r1, r2 - 1
392 tmp = work( inds+i ) + work( indp+i )
393 IF( tmp.EQ.zero )
394 $ tmp = eps*work( inds+i )
395 IF( abs( tmp ).LE.abs( mingma ) ) THEN
396 mingma = tmp
397 r = i + 1
398 END IF
399 110 CONTINUE
400*
401* Compute the FP vector: solve N^T v = e_r
402*
403 isuppz( 1 ) = b1
404 isuppz( 2 ) = bn
405 z( r ) = cone
406 ztz = one
407*
408* Compute the FP vector upwards from R
409*
410 IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
411 DO 210 i = r-1, b1, -1
412 z( i ) = -( work( indlpl+i )*z( i+1 ) )
413 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
414 $ THEN
415 z( i ) = zero
416 isuppz( 1 ) = i + 1
417 GOTO 220
418 ENDIF
419 ztz = ztz + real( z( i )*z( i ) )
420 210 CONTINUE
421 220 CONTINUE
422 ELSE
423* Run slower loop if NaN occurred.
424 DO 230 i = r - 1, b1, -1
425 IF( z( i+1 ).EQ.zero ) THEN
426 z( i ) = -( ld( i+1 ) / ld( i ) )*z( i+2 )
427 ELSE
428 z( i ) = -( work( indlpl+i )*z( i+1 ) )
429 END IF
430 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
431 $ THEN
432 z( i ) = zero
433 isuppz( 1 ) = i + 1
434 GO TO 240
435 END IF
436 ztz = ztz + real( z( i )*z( i ) )
437 230 CONTINUE
438 240 CONTINUE
439 ENDIF
440
441* Compute the FP vector downwards from R in blocks of size BLKSIZ
442 IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
443 DO 250 i = r, bn-1
444 z( i+1 ) = -( work( indumn+i )*z( i ) )
445 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
446 $ THEN
447 z( i+1 ) = zero
448 isuppz( 2 ) = i
449 GO TO 260
450 END IF
451 ztz = ztz + real( z( i+1 )*z( i+1 ) )
452 250 CONTINUE
453 260 CONTINUE
454 ELSE
455* Run slower loop if NaN occurred.
456 DO 270 i = r, bn - 1
457 IF( z( i ).EQ.zero ) THEN
458 z( i+1 ) = -( ld( i-1 ) / ld( i ) )*z( i-1 )
459 ELSE
460 z( i+1 ) = -( work( indumn+i )*z( i ) )
461 END IF
462 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
463 $ THEN
464 z( i+1 ) = zero
465 isuppz( 2 ) = i
466 GO TO 280
467 END IF
468 ztz = ztz + real( z( i+1 )*z( i+1 ) )
469 270 CONTINUE
470 280 CONTINUE
471 END IF
472*
473* Compute quantities for convergence test
474*
475 tmp = one / ztz
476 nrminv = sqrt( tmp )
477 resid = abs( mingma )*nrminv
478 rqcorr = mingma*tmp
479*
480*
481 RETURN
482*
483* End of CLAR1V
484*
485 END
subroutine clar1v(n, b1, bn, lambda, d, l, ld, lld, pivmin, gaptol, z, wantnc, negcnt, ztz, mingma, r, isuppz, nrminv, resid, rqcorr, work)
CLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition clar1v.f:230