LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zlar1v ( 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,
complex*16, 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 
)

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

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

Purpose: 
 ZLAR1V 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 COMPLEX*16 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.
Date
September 2012
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 232 of file zlar1v.f.

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

Here is the caller graph for this function: