floraachy
/
OpenBLAS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
OpenBLAS/lapack-netlib/SRC/dlar1v.f

487 lines
14 KiB
Fortran
Raw Blame History

*> \brief \b DLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAR1V + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlar1v.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlar1v.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlar1v.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
* PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
* R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
*
* .. Scalar Arguments ..
* LOGICAL WANTNC
* INTEGER B1, BN, N, NEGCNT, R
* DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
* $ RQCORR, ZTZ
* ..
* .. Array Arguments ..
* INTEGER ISUPPZ( * )
* DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
* $ WORK( * )
* DOUBLE PRECISION Z( * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The order of the matrix L D L**T.
*> \endverbatim
*>
*> \param[in] B1
*> \verbatim
*> B1 is INTEGER
*> First index of the submatrix of L D L**T.
*> \endverbatim
*>
*> \param[in] BN
*> \verbatim
*> BN is INTEGER
*> Last index of the submatrix of L D L**T.
*> \endverbatim
*>
*> \param[in] LAMBDA
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in] L
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in] D
*> \verbatim
*> D is DOUBLE PRECISION array, dimension (N)
*> The n diagonal elements of the diagonal matrix D.
*> \endverbatim
*>
*> \param[in] LD
*> \verbatim
*> LD is DOUBLE PRECISION array, dimension (N-1)
*> The n-1 elements L(i)*D(i).
*> \endverbatim
*>
*> \param[in] LLD
*> \verbatim
*> LLD is DOUBLE PRECISION array, dimension (N-1)
*> The n-1 elements L(i)*L(i)*D(i).
*> \endverbatim
*>
*> \param[in] PIVMIN
*> \verbatim
*> PIVMIN is DOUBLE PRECISION
*> The minimum pivot in the Sturm sequence.
*> \endverbatim
*>
*> \param[in] GAPTOL
*> \verbatim
*> GAPTOL is DOUBLE PRECISION
*> Tolerance that indicates when eigenvector entries are negligible
*> w.r.t. their contribution to the residual.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[in] WANTNC
*> \verbatim
*> WANTNC is LOGICAL
*> Specifies whether NEGCNT has to be computed.
*> \endverbatim
*>
*> \param[out] NEGCNT
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[out] ZTZ
*> \verbatim
*> ZTZ is DOUBLE PRECISION
*> The square of the 2-norm of Z.
*> \endverbatim
*>
*> \param[out] MINGMA
*> \verbatim
*> MINGMA is DOUBLE PRECISION
*> The reciprocal of the largest (in magnitude) diagonal
*> element of the inverse of L D L**T - sigma I.
*> \endverbatim
*>
*> \param[in,out] R
*> \verbatim
*> 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.
*> \endverbatim
*>
*> \param[out] ISUPPZ
*> \verbatim
*> 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 ).
*> \endverbatim
*>
*> \param[out] NRMINV
*> \verbatim
*> NRMINV is DOUBLE PRECISION
*> NRMINV = 1/SQRT( ZTZ )
*> \endverbatim
*>
*> \param[out] RESID
*> \verbatim
*> RESID is DOUBLE PRECISION
*> The residual of the FP vector.
*> RESID = ABS( MINGMA )/SQRT( ZTZ )
*> \endverbatim
*>
*> \param[out] RQCORR
*> \verbatim
*> RQCORR is DOUBLE PRECISION
*> The Rayleigh Quotient correction to LAMBDA.
*> RQCORR = MINGMA*TMP
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*> WORK is DOUBLE PRECISION array, dimension (4*N)
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date December 2016
*
*> \ingroup doubleOTHERauxiliary
*
*> \par Contributors:
* ==================
*>
*> Beresford Parlett, University of California, Berkeley, USA \n
*> Jim Demmel, University of California, Berkeley, USA \n
*> Inderjit Dhillon, University of Texas, Austin, USA \n
*> Osni Marques, LBNL/NERSC, USA \n
*> Christof Voemel, University of California, Berkeley, USA
*
* =====================================================================
SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
$ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
$ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
*
* -- LAPACK auxiliary routine (version 3.7.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
* December 2016
*
* .. Scalar Arguments ..
LOGICAL WANTNC
INTEGER B1, BN, N, NEGCNT, R
DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
$ RQCORR, ZTZ
* ..
* .. Array Arguments ..
INTEGER ISUPPZ( * )
DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
$ WORK( * )
DOUBLE PRECISION Z( * )
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
* ..
* .. Local Scalars ..
LOGICAL SAWNAN1, SAWNAN2
INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
$ R2
DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
* ..
* .. External Functions ..
LOGICAL DISNAN
DOUBLE PRECISION DLAMCH
EXTERNAL DISNAN, DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS
* ..
* .. Executable Statements ..
*
EPS = DLAMCH( 'Precision' )
IF( R.EQ.0 ) THEN
R1 = B1
R2 = BN
ELSE
R1 = R
R2 = R
END IF
* Storage for LPLUS
INDLPL = 0
* Storage for UMINUS
INDUMN = N
INDS = 2*N + 1
INDP = 3*N + 1
IF( B1.EQ.1 ) THEN
WORK( INDS ) = ZERO
ELSE
WORK( INDS+B1-1 ) = LLD( B1-1 )
END IF
*
* Compute the stationary transform (using the differential form)
* until the index R2.
*
SAWNAN1 = .FALSE.
NEG1 = 0
S = WORK( INDS+B1-1 ) - LAMBDA
DO 50 I = B1, R1 - 1
DPLUS = D( I ) + S
WORK( INDLPL+I ) = LD( I ) / DPLUS
IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
S = WORK( INDS+I ) - LAMBDA
50 CONTINUE
SAWNAN1 = DISNAN( S )
IF( SAWNAN1 ) GOTO 60
DO 51 I = R1, R2 - 1
DPLUS = D( I ) + S
WORK( INDLPL+I ) = LD( I ) / DPLUS
WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
S = WORK( INDS+I ) - LAMBDA
51 CONTINUE
SAWNAN1 = DISNAN( S )
*
60 CONTINUE
IF( SAWNAN1 ) THEN
* Runs a slower version of the above loop if a NaN is detected
NEG1 = 0
S = WORK( INDS+B1-1 ) - LAMBDA
DO 70 I = B1, R1 - 1
DPLUS = D( I ) + S
IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
WORK( INDLPL+I ) = LD( I ) / DPLUS
IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
IF( WORK( INDLPL+I ).EQ.ZERO )
$ WORK( INDS+I ) = LLD( I )
S = WORK( INDS+I ) - LAMBDA
70 CONTINUE
DO 71 I = R1, R2 - 1
DPLUS = D( I ) + S
IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
WORK( INDLPL+I ) = LD( I ) / DPLUS
WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
IF( WORK( INDLPL+I ).EQ.ZERO )
$ WORK( INDS+I ) = LLD( I )
S = WORK( INDS+I ) - LAMBDA
71 CONTINUE
END IF
*
* Compute the progressive transform (using the differential form)
* until the index R1
*
SAWNAN2 = .FALSE.
NEG2 = 0
WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
DO 80 I = BN - 1, R1, -1
DMINUS = LLD( I ) + WORK( INDP+I )
TMP = D( I ) / DMINUS
IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
WORK( INDUMN+I ) = L( I )*TMP
WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
80 CONTINUE
TMP = WORK( INDP+R1-1 )
SAWNAN2 = DISNAN( TMP )
IF( SAWNAN2 ) THEN
* Runs a slower version of the above loop if a NaN is detected
NEG2 = 0
DO 100 I = BN-1, R1, -1
DMINUS = LLD( I ) + WORK( INDP+I )
IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
TMP = D( I ) / DMINUS
IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
WORK( INDUMN+I ) = L( I )*TMP
WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
IF( TMP.EQ.ZERO )
$ WORK( INDP+I-1 ) = D( I ) - LAMBDA
100 CONTINUE
END IF
*
* Find the index (from R1 to R2) of the largest (in magnitude)
* diagonal element of the inverse
*
MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
IF( WANTNC ) THEN
NEGCNT = NEG1 + NEG2
ELSE
NEGCNT = -1
ENDIF
IF( ABS(MINGMA).EQ.ZERO )
$ MINGMA = EPS*WORK( INDS+R1-1 )
R = R1
DO 110 I = R1, R2 - 1
TMP = WORK( INDS+I ) + WORK( INDP+I )
IF( TMP.EQ.ZERO )
$ TMP = EPS*WORK( INDS+I )
IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
MINGMA = TMP
R = I + 1
END IF
110 CONTINUE
*
* Compute the FP vector: solve N^T v = e_r
*
ISUPPZ( 1 ) = B1
ISUPPZ( 2 ) = BN
Z( R ) = ONE
ZTZ = ONE
*
* Compute the FP vector upwards from R
*
IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
DO 210 I = R-1, B1, -1
Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
$ THEN
Z( I ) = ZERO
ISUPPZ( 1 ) = I + 1
GOTO 220
ENDIF
ZTZ = ZTZ + Z( I )*Z( I )
210 CONTINUE
220 CONTINUE
ELSE
* Run slower loop if NaN occurred.
DO 230 I = R - 1, B1, -1
IF( Z( I+1 ).EQ.ZERO ) THEN
Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
ELSE
Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
END IF
IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
$ THEN
Z( I ) = ZERO
ISUPPZ( 1 ) = I + 1
GO TO 240
END IF
ZTZ = ZTZ + Z( I )*Z( I )
230 CONTINUE
240 CONTINUE
ENDIF
* Compute the FP vector downwards from R in blocks of size BLKSIZ
IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
DO 250 I = R, BN-1
Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
$ THEN
Z( I+1 ) = ZERO
ISUPPZ( 2 ) = I
GO TO 260
END IF
ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
250 CONTINUE
260 CONTINUE
ELSE
* Run slower loop if NaN occurred.
DO 270 I = R, BN - 1
IF( Z( I ).EQ.ZERO ) THEN
Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
ELSE
Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
END IF
IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
$ THEN
Z( I+1 ) = ZERO
ISUPPZ( 2 ) = I
GO TO 280
END IF
ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
270 CONTINUE
280 CONTINUE
END IF
*
* Compute quantities for convergence test
*
TMP = ONE / ZTZ
NRMINV = SQRT( TMP )
RESID = ABS( MINGMA )*NRMINV
RQCORR = MINGMA*TMP
*
*
RETURN
*
* End of DLAR1V
*
END
Reference in New Issue View Git Blame Copy Permalink