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

◆ clar1v()

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

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

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

Purpose:
!>
!> CLAR1V 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 COMPLEX 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 225 of file clar1v.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 REAL GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
237 $ RQCORR, ZTZ
238* ..
239* .. Array Arguments ..
240 INTEGER ISUPPZ( * )
241 REAL D( * ), L( * ), LD( * ), LLD( * ),
242 $ WORK( * )
243 COMPLEX Z( * )
244* ..
245*
246* =====================================================================
247*
248* .. Parameters ..
249 REAL ZERO, ONE
250 parameter( zero = 0.0e0, one = 1.0e0 )
251 COMPLEX CONE
252 parameter( cone = ( 1.0e0, 0.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, real
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 ) = cone
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 + real( 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 + real( 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 + real( 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 + real( 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 CLAR1V
482*
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:57
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
Here is the caller graph for this function: