 LAPACK  3.8.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.

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)`
Date
December 2016
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 slar1v.f.

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