LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlar1v()

subroutine dlar1v ( integer n,
integer b1,
integer bn,
double precision lambda,
double precision, dimension( * ) d,
double precision, dimension( * ) l,
double precision, dimension( * ) ld,
double precision, dimension( * ) lld,
double precision pivmin,
double precision gaptol,
double precision, dimension( * ) z,
logical wantnc,
integer negcnt,
double precision ztz,
double precision mingma,
integer r,
integer, dimension( * ) isuppz,
double precision nrminv,
double precision resid,
double precision rqcorr,
double precision, dimension( * ) work )

DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.

Download DLAR1V + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLAR1V computes the (scaled) r-th column of the inverse of
!> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
!> L D L**T - sigma I. When sigma is close to an eigenvalue, the
!> computed vector is an accurate eigenvector. Usually, r corresponds
!> to the index where the eigenvector is largest in magnitude.
!> The following steps accomplish this computation :
!> (a) Stationary qd transform,  L D L**T - sigma I = L(+) D(+) L(+)**T,
!> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
!> (c) Computation of the diagonal elements of the inverse of
!>     L D L**T - sigma I by combining the above transforms, and choosing
!>     r as the index where the diagonal of the inverse is (one of the)
!>     largest in magnitude.
!> (d) Computation of the (scaled) r-th column of the inverse using the
!>     twisted factorization obtained by combining the top part of the
!>     the stationary and the bottom part of the progressive transform.
!>
Parameters
[in]N
!>          N is INTEGER
!>           The order of the matrix L D L**T.
!>
[in]B1
!>          B1 is INTEGER
!>           First index of the submatrix of L D L**T.
!>
[in]BN
!>          BN is INTEGER
!>           Last index of the submatrix of L D L**T.
!>
[in]LAMBDA
!>          LAMBDA is DOUBLE PRECISION
!>           The shift. In order to compute an accurate eigenvector,
!>           LAMBDA should be a good approximation to an eigenvalue
!>           of L D L**T.
!>
[in]L
!>          L is DOUBLE PRECISION array, dimension (N-1)
!>           The (n-1) subdiagonal elements of the unit bidiagonal matrix
!>           L, in elements 1 to N-1.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>           The n diagonal elements of the diagonal matrix D.
!>
[in]LD
!>          LD is DOUBLE PRECISION array, dimension (N-1)
!>           The n-1 elements L(i)*D(i).
!>
[in]LLD
!>          LLD is DOUBLE PRECISION array, dimension (N-1)
!>           The n-1 elements L(i)*L(i)*D(i).
!>
[in]PIVMIN
!>          PIVMIN is DOUBLE PRECISION
!>           The minimum pivot in the Sturm sequence.
!>
[in]GAPTOL
!>          GAPTOL is DOUBLE PRECISION
!>           Tolerance that indicates when eigenvector entries are negligible
!>           w.r.t. their contribution to the residual.
!>
[in,out]Z
!>          Z is DOUBLE PRECISION array, dimension (N)
!>           On input, all entries of Z must be set to 0.
!>           On output, Z contains the (scaled) r-th column of the
!>           inverse. The scaling is such that Z(R) equals 1.
!>
[in]WANTNC
!>          WANTNC is LOGICAL
!>           Specifies whether NEGCNT has to be computed.
!>
[out]NEGCNT
!>          NEGCNT is INTEGER
!>           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
!>           in the  matrix factorization L D L**T, and NEGCNT = -1 otherwise.
!>
[out]ZTZ
!>          ZTZ is DOUBLE PRECISION
!>           The square of the 2-norm of Z.
!>
[out]MINGMA
!>          MINGMA is DOUBLE PRECISION
!>           The reciprocal of the largest (in magnitude) diagonal
!>           element of the inverse of L D L**T - sigma I.
!>
[in,out]R
!>          R is INTEGER
!>           The twist index for the twisted factorization used to
!>           compute Z.
!>           On input, 0 <= R <= N. If R is input as 0, R is set to
!>           the index where (L D L**T - sigma I)^{-1} is largest
!>           in magnitude. If 1 <= R <= N, R is unchanged.
!>           On output, R contains the twist index used to compute Z.
!>           Ideally, R designates the position of the maximum entry in the
!>           eigenvector.
!>
[out]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension (2)
!>           The support of the vector in Z, i.e., the vector Z is
!>           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
!>
[out]NRMINV
!>          NRMINV is DOUBLE PRECISION
!>           NRMINV = 1/SQRT( ZTZ )
!>
[out]RESID
!>          RESID is DOUBLE PRECISION
!>           The residual of the FP vector.
!>           RESID = ABS( MINGMA )/SQRT( ZTZ )
!>
[out]RQCORR
!>          RQCORR is DOUBLE PRECISION
!>           The Rayleigh Quotient correction to LAMBDA.
!>           RQCORR = MINGMA*TMP
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (4*N)
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 225 of file dlar1v.f.

228*
229* -- LAPACK auxiliary routine --
230* -- LAPACK is a software package provided by Univ. of Tennessee, --
231* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232*
233* .. Scalar Arguments ..
234 LOGICAL WANTNC
235 INTEGER B1, BN, N, NEGCNT, R
236 DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
237 $ RQCORR, ZTZ
238* ..
239* .. Array Arguments ..
240 INTEGER ISUPPZ( * )
241 DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
242 $ WORK( * )
243 DOUBLE PRECISION Z( * )
244* ..
245*
246* =====================================================================
247*
248* .. Parameters ..
249 DOUBLE PRECISION ZERO, ONE
250 parameter( zero = 0.0d0, one = 1.0d0 )
251
252* ..
253* .. Local Scalars ..
254 LOGICAL SAWNAN1, SAWNAN2
255 INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
256 $ R2
257 DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
258* ..
259* .. External Functions ..
260 LOGICAL DISNAN
261 DOUBLE PRECISION DLAMCH
262 EXTERNAL disnan, dlamch
263* ..
264* .. Intrinsic Functions ..
265 INTRINSIC abs
266* ..
267* .. Executable Statements ..
268*
269 eps = dlamch( 'Precision' )
270
271
272 IF( r.EQ.0 ) THEN
273 r1 = b1
274 r2 = bn
275 ELSE
276 r1 = r
277 r2 = r
278 END IF
279
280* Storage for LPLUS
281 indlpl = 0
282* Storage for UMINUS
283 indumn = n
284 inds = 2*n + 1
285 indp = 3*n + 1
286
287 IF( b1.EQ.1 ) THEN
288 work( inds ) = zero
289 ELSE
290 work( inds+b1-1 ) = lld( b1-1 )
291 END IF
292
293*
294* Compute the stationary transform (using the differential form)
295* until the index R2.
296*
297 sawnan1 = .false.
298 neg1 = 0
299 s = work( inds+b1-1 ) - lambda
300 DO 50 i = b1, r1 - 1
301 dplus = d( i ) + s
302 work( indlpl+i ) = ld( i ) / dplus
303 IF(dplus.LT.zero) neg1 = neg1 + 1
304 work( inds+i ) = s*work( indlpl+i )*l( i )
305 s = work( inds+i ) - lambda
306 50 CONTINUE
307 sawnan1 = disnan( s )
308 IF( sawnan1 ) GOTO 60
309 DO 51 i = r1, r2 - 1
310 dplus = d( i ) + s
311 work( indlpl+i ) = ld( i ) / dplus
312 work( inds+i ) = s*work( indlpl+i )*l( i )
313 s = work( inds+i ) - lambda
314 51 CONTINUE
315 sawnan1 = disnan( s )
316*
317 60 CONTINUE
318 IF( sawnan1 ) THEN
319* Runs a slower version of the above loop if a NaN is detected
320 neg1 = 0
321 s = work( inds+b1-1 ) - lambda
322 DO 70 i = b1, r1 - 1
323 dplus = d( i ) + s
324 IF(abs(dplus).LT.pivmin) dplus = -pivmin
325 work( indlpl+i ) = ld( i ) / dplus
326 IF(dplus.LT.zero) neg1 = neg1 + 1
327 work( inds+i ) = s*work( indlpl+i )*l( i )
328 IF( work( indlpl+i ).EQ.zero )
329 $ work( inds+i ) = lld( i )
330 s = work( inds+i ) - lambda
331 70 CONTINUE
332 DO 71 i = r1, r2 - 1
333 dplus = d( i ) + s
334 IF(abs(dplus).LT.pivmin) dplus = -pivmin
335 work( indlpl+i ) = ld( i ) / dplus
336 work( inds+i ) = s*work( indlpl+i )*l( i )
337 IF( work( indlpl+i ).EQ.zero )
338 $ work( inds+i ) = lld( i )
339 s = work( inds+i ) - lambda
340 71 CONTINUE
341 END IF
342*
343* Compute the progressive transform (using the differential form)
344* until the index R1
345*
346 sawnan2 = .false.
347 neg2 = 0
348 work( indp+bn-1 ) = d( bn ) - lambda
349 DO 80 i = bn - 1, r1, -1
350 dminus = lld( i ) + work( indp+i )
351 tmp = d( i ) / dminus
352 IF(dminus.LT.zero) neg2 = neg2 + 1
353 work( indumn+i ) = l( i )*tmp
354 work( indp+i-1 ) = work( indp+i )*tmp - lambda
355 80 CONTINUE
356 tmp = work( indp+r1-1 )
357 sawnan2 = disnan( tmp )
358
359 IF( sawnan2 ) THEN
360* Runs a slower version of the above loop if a NaN is detected
361 neg2 = 0
362 DO 100 i = bn-1, r1, -1
363 dminus = lld( i ) + work( indp+i )
364 IF(abs(dminus).LT.pivmin) dminus = -pivmin
365 tmp = d( i ) / dminus
366 IF(dminus.LT.zero) neg2 = neg2 + 1
367 work( indumn+i ) = l( i )*tmp
368 work( indp+i-1 ) = work( indp+i )*tmp - lambda
369 IF( tmp.EQ.zero )
370 $ work( indp+i-1 ) = d( i ) - lambda
371 100 CONTINUE
372 END IF
373*
374* Find the index (from R1 to R2) of the largest (in magnitude)
375* diagonal element of the inverse
376*
377 mingma = work( inds+r1-1 ) + work( indp+r1-1 )
378 IF( mingma.LT.zero ) neg1 = neg1 + 1
379 IF( wantnc ) THEN
380 negcnt = neg1 + neg2
381 ELSE
382 negcnt = -1
383 ENDIF
384 IF( abs(mingma).EQ.zero )
385 $ mingma = eps*work( inds+r1-1 )
386 r = r1
387 DO 110 i = r1, r2 - 1
388 tmp = work( inds+i ) + work( indp+i )
389 IF( tmp.EQ.zero )
390 $ tmp = eps*work( inds+i )
391 IF( abs( tmp ).LE.abs( mingma ) ) THEN
392 mingma = tmp
393 r = i + 1
394 END IF
395 110 CONTINUE
396*
397* Compute the FP vector: solve N^T v = e_r
398*
399 isuppz( 1 ) = b1
400 isuppz( 2 ) = bn
401 z( r ) = one
402 ztz = one
403*
404* Compute the FP vector upwards from R
405*
406 IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
407 DO 210 i = r-1, b1, -1
408 z( i ) = -( work( indlpl+i )*z( i+1 ) )
409 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
410 $ THEN
411 z( i ) = zero
412 isuppz( 1 ) = i + 1
413 GOTO 220
414 ENDIF
415 ztz = ztz + z( i )*z( i )
416 210 CONTINUE
417 220 CONTINUE
418 ELSE
419* Run slower loop if NaN occurred.
420 DO 230 i = r - 1, b1, -1
421 IF( z( i+1 ).EQ.zero ) THEN
422 z( i ) = -( ld( i+1 ) / ld( i ) )*z( i+2 )
423 ELSE
424 z( i ) = -( work( indlpl+i )*z( i+1 ) )
425 END IF
426 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
427 $ THEN
428 z( i ) = zero
429 isuppz( 1 ) = i + 1
430 GO TO 240
431 END IF
432 ztz = ztz + z( i )*z( i )
433 230 CONTINUE
434 240 CONTINUE
435 ENDIF
436
437* Compute the FP vector downwards from R in blocks of size BLKSIZ
438 IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
439 DO 250 i = r, bn-1
440 z( i+1 ) = -( work( indumn+i )*z( i ) )
441 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
442 $ THEN
443 z( i+1 ) = zero
444 isuppz( 2 ) = i
445 GO TO 260
446 END IF
447 ztz = ztz + z( i+1 )*z( i+1 )
448 250 CONTINUE
449 260 CONTINUE
450 ELSE
451* Run slower loop if NaN occurred.
452 DO 270 i = r, bn - 1
453 IF( z( i ).EQ.zero ) THEN
454 z( i+1 ) = -( ld( i-1 ) / ld( i ) )*z( i-1 )
455 ELSE
456 z( i+1 ) = -( work( indumn+i )*z( i ) )
457 END IF
458 IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
459 $ THEN
460 z( i+1 ) = zero
461 isuppz( 2 ) = i
462 GO TO 280
463 END IF
464 ztz = ztz + z( i+1 )*z( i+1 )
465 270 CONTINUE
466 280 CONTINUE
467 END IF
468*
469* Compute quantities for convergence test
470*
471 tmp = one / ztz
472 nrminv = sqrt( tmp )
473 resid = abs( mingma )*nrminv
474 rqcorr = mingma*tmp
475*
476*
477 RETURN
478*
479* End of DLAR1V
480*
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:57
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
Here is the caller graph for this function: