LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slar1v()

subroutine slar1v ( integer  N,
integer  B1,
integer  BN,
real  LAMBDA,
real, dimension( * )  D,
real, dimension( * )  L,
real, dimension( * )  LD,
real, dimension( * )  LLD,
real  PIVMIN,
real  GAPTOL,
real, dimension( * )  Z,
logical  WANTNC,
integer  NEGCNT,
real  ZTZ,
real  MINGMA,
integer  R,
integer, dimension( * )  ISUPPZ,
real  NRMINV,
real  RESID,
real  RQCORR,
real, dimension( * )  WORK 
)

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

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

Purpose:
 SLAR1V 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 REAL
           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 REAL 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 REAL array, dimension (N)
           The n diagonal elements of the diagonal matrix D.
[in]LD
          LD is REAL array, dimension (N-1)
           The n-1 elements L(i)*D(i).
[in]LLD
          LLD is REAL array, dimension (N-1)
           The n-1 elements L(i)*L(i)*D(i).
[in]PIVMIN
          PIVMIN is REAL
           The minimum pivot in the Sturm sequence.
[in]GAPTOL
          GAPTOL is REAL
           Tolerance that indicates when eigenvector entries are negligible
           w.r.t. their contribution to the residual.
[in,out]Z
          Z is REAL 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 REAL
           The square of the 2-norm of Z.
[out]MINGMA
          MINGMA is REAL
           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 REAL
           NRMINV = 1/SQRT( ZTZ )
[out]RESID
          RESID is REAL
           The residual of the FP vector.
           RESID = ABS( MINGMA )/SQRT( ZTZ )
[out]RQCORR
          RQCORR is REAL
           The Rayleigh Quotient correction to LAMBDA.
           RQCORR = MINGMA*TMP
[out]WORK
          WORK is REAL 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 227 of file slar1v.f.

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