LAPACK: SRC/dlar1v.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
00002      $           PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
00003      $           R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
00004 *
00005 *  -- LAPACK auxiliary routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       LOGICAL            WANTNC
00012       INTEGER   B1, BN, N, NEGCNT, R
00013       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
00014      $                   RQCORR, ZTZ
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            ISUPPZ( * )
00018       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
00019      $                  WORK( * )
00020       DOUBLE PRECISION Z( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  DLAR1V computes the (scaled) r-th column of the inverse of
00027 *  the sumbmatrix in rows B1 through BN of the tridiagonal matrix
00028 *  L D L^T - sigma I. When sigma is close to an eigenvalue, the
00029 *  computed vector is an accurate eigenvector. Usually, r corresponds
00030 *  to the index where the eigenvector is largest in magnitude.
00031 *  The following steps accomplish this computation :
00032 *  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T,
00033 *  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
00034 *  (c) Computation of the diagonal elements of the inverse of
00035 *      L D L^T - sigma I by combining the above transforms, and choosing
00036 *      r as the index where the diagonal of the inverse is (one of the)
00037 *      largest in magnitude.
00038 *  (d) Computation of the (scaled) r-th column of the inverse using the
00039 *      twisted factorization obtained by combining the top part of the
00040 *      the stationary and the bottom part of the progressive transform.
00041 *
00042 *  Arguments
00043 *  =========
00044 *
00045 *  N        (input) INTEGER
00046 *           The order of the matrix L D L^T.
00047 *
00048 *  B1       (input) INTEGER
00049 *           First index of the submatrix of L D L^T.
00050 *
00051 *  BN       (input) INTEGER
00052 *           Last index of the submatrix of L D L^T.
00053 *
00054 *  LAMBDA    (input) DOUBLE PRECISION
00055 *           The shift. In order to compute an accurate eigenvector,
00056 *           LAMBDA should be a good approximation to an eigenvalue
00057 *           of L D L^T.
00058 *
00059 *  L        (input) DOUBLE PRECISION array, dimension (N-1)
00060 *           The (n-1) subdiagonal elements of the unit bidiagonal matrix
00061 *           L, in elements 1 to N-1.
00062 *
00063 *  D        (input) DOUBLE PRECISION array, dimension (N)
00064 *           The n diagonal elements of the diagonal matrix D.
00065 *
00066 *  LD       (input) DOUBLE PRECISION array, dimension (N-1)
00067 *           The n-1 elements L(i)*D(i).
00068 *
00069 *  LLD      (input) DOUBLE PRECISION array, dimension (N-1)
00070 *           The n-1 elements L(i)*L(i)*D(i).
00071 *
00072 *  PIVMIN   (input) DOUBLE PRECISION
00073 *           The minimum pivot in the Sturm sequence.
00074 *
00075 *  GAPTOL   (input) DOUBLE PRECISION
00076 *           Tolerance that indicates when eigenvector entries are negligible
00077 *           w.r.t. their contribution to the residual.
00078 *
00079 *  Z        (input/output) DOUBLE PRECISION array, dimension (N)
00080 *           On input, all entries of Z must be set to 0.
00081 *           On output, Z contains the (scaled) r-th column of the
00082 *           inverse. The scaling is such that Z(R) equals 1.
00083 *
00084 *  WANTNC   (input) LOGICAL
00085 *           Specifies whether NEGCNT has to be computed.
00086 *
00087 *  NEGCNT   (output) INTEGER
00088 *           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
00089 *           in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise.
00090 *
00091 *  ZTZ      (output) DOUBLE PRECISION
00092 *           The square of the 2-norm of Z.
00093 *
00094 *  MINGMA   (output) DOUBLE PRECISION
00095 *           The reciprocal of the largest (in magnitude) diagonal
00096 *           element of the inverse of L D L^T - sigma I.
00097 *
00098 *  R        (input/output) INTEGER
00099 *           The twist index for the twisted factorization used to
00100 *           compute Z.
00101 *           On input, 0 <= R <= N. If R is input as 0, R is set to
00102 *           the index where (L D L^T - sigma I)^{-1} is largest
00103 *           in magnitude. If 1 <= R <= N, R is unchanged.
00104 *           On output, R contains the twist index used to compute Z.
00105 *           Ideally, R designates the position of the maximum entry in the
00106 *           eigenvector.
00107 *
00108 *  ISUPPZ   (output) INTEGER array, dimension (2)
00109 *           The support of the vector in Z, i.e., the vector Z is
00110 *           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
00111 *
00112 *  NRMINV   (output) DOUBLE PRECISION
00113 *           NRMINV = 1/SQRT( ZTZ )
00114 *
00115 *  RESID    (output) DOUBLE PRECISION
00116 *           The residual of the FP vector.
00117 *           RESID = ABS( MINGMA )/SQRT( ZTZ )
00118 *
00119 *  RQCORR   (output) DOUBLE PRECISION
00120 *           The Rayleigh Quotient correction to LAMBDA.
00121 *           RQCORR = MINGMA*TMP
00122 *
00123 *  WORK     (workspace) DOUBLE PRECISION array, dimension (4*N)
00124 *
00125 *  Further Details
00126 *  ===============
00127 *
00128 *  Based on contributions by
00129 *     Beresford Parlett, University of California, Berkeley, USA
00130 *     Jim Demmel, University of California, Berkeley, USA
00131 *     Inderjit Dhillon, University of Texas, Austin, USA
00132 *     Osni Marques, LBNL/NERSC, USA
00133 *     Christof Voemel, University of California, Berkeley, USA
00134 *
00135 *  =====================================================================
00136 *
00137 *     .. Parameters ..
00138       DOUBLE PRECISION   ZERO, ONE
00139       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00140 
00141 *     ..
00142 *     .. Local Scalars ..
00143       LOGICAL            SAWNAN1, SAWNAN2
00144       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
00145      $                   R2
00146       DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
00147 *     ..
00148 *     .. External Functions ..
00149       LOGICAL DISNAN
00150       DOUBLE PRECISION   DLAMCH
00151       EXTERNAL           DISNAN, DLAMCH
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC          ABS
00155 *     ..
00156 *     .. Executable Statements ..
00157 *
00158       EPS = DLAMCH( 'Precision' )
00159 
00160 
00161       IF( R.EQ.0 ) THEN
00162          R1 = B1
00163          R2 = BN
00164       ELSE
00165          R1 = R
00166          R2 = R
00167       END IF
00168 
00169 *     Storage for LPLUS
00170       INDLPL = 0
00171 *     Storage for UMINUS
00172       INDUMN = N
00173       INDS = 2*N + 1
00174       INDP = 3*N + 1
00175 
00176       IF( B1.EQ.1 ) THEN
00177          WORK( INDS ) = ZERO
00178       ELSE
00179          WORK( INDS+B1-1 ) = LLD( B1-1 )
00180       END IF
00181 
00182 *
00183 *     Compute the stationary transform (using the differential form)
00184 *     until the index R2.
00185 *
00186       SAWNAN1 = .FALSE.
00187       NEG1 = 0
00188       S = WORK( INDS+B1-1 ) - LAMBDA
00189       DO 50 I = B1, R1 - 1
00190          DPLUS = D( I ) + S
00191          WORK( INDLPL+I ) = LD( I ) / DPLUS
00192          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00193          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00194          S = WORK( INDS+I ) - LAMBDA
00195  50   CONTINUE
00196       SAWNAN1 = DISNAN( S )
00197       IF( SAWNAN1 ) GOTO 60
00198       DO 51 I = R1, R2 - 1
00199          DPLUS = D( I ) + S
00200          WORK( INDLPL+I ) = LD( I ) / DPLUS
00201          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00202          S = WORK( INDS+I ) - LAMBDA
00203  51   CONTINUE
00204       SAWNAN1 = DISNAN( S )
00205 *
00206  60   CONTINUE
00207       IF( SAWNAN1 ) THEN
00208 *        Runs a slower version of the above loop if a NaN is detected
00209          NEG1 = 0
00210          S = WORK( INDS+B1-1 ) - LAMBDA
00211          DO 70 I = B1, R1 - 1
00212             DPLUS = D( I ) + S
00213             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00214             WORK( INDLPL+I ) = LD( I ) / DPLUS
00215             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00216             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00217             IF( WORK( INDLPL+I ).EQ.ZERO )
00218      $                      WORK( INDS+I ) = LLD( I )
00219             S = WORK( INDS+I ) - LAMBDA
00220  70      CONTINUE
00221          DO 71 I = R1, R2 - 1
00222             DPLUS = D( I ) + S
00223             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00224             WORK( INDLPL+I ) = LD( I ) / DPLUS
00225             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00226             IF( WORK( INDLPL+I ).EQ.ZERO )
00227      $                      WORK( INDS+I ) = LLD( I )
00228             S = WORK( INDS+I ) - LAMBDA
00229  71      CONTINUE
00230       END IF
00231 *
00232 *     Compute the progressive transform (using the differential form)
00233 *     until the index R1
00234 *
00235       SAWNAN2 = .FALSE.
00236       NEG2 = 0
00237       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
00238       DO 80 I = BN - 1, R1, -1
00239          DMINUS = LLD( I ) + WORK( INDP+I )
00240          TMP = D( I ) / DMINUS
00241          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00242          WORK( INDUMN+I ) = L( I )*TMP
00243          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00244  80   CONTINUE
00245       TMP = WORK( INDP+R1-1 )
00246       SAWNAN2 = DISNAN( TMP )
00247 
00248       IF( SAWNAN2 ) THEN
00249 *        Runs a slower version of the above loop if a NaN is detected
00250          NEG2 = 0
00251          DO 100 I = BN-1, R1, -1
00252             DMINUS = LLD( I ) + WORK( INDP+I )
00253             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
00254             TMP = D( I ) / DMINUS
00255             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00256             WORK( INDUMN+I ) = L( I )*TMP
00257             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00258             IF( TMP.EQ.ZERO )
00259      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
00260  100     CONTINUE
00261       END IF
00262 *
00263 *     Find the index (from R1 to R2) of the largest (in magnitude)
00264 *     diagonal element of the inverse
00265 *
00266       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
00267       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
00268       IF( WANTNC ) THEN
00269          NEGCNT = NEG1 + NEG2
00270       ELSE
00271          NEGCNT = -1
00272       ENDIF
00273       IF( ABS(MINGMA).EQ.ZERO )
00274      $   MINGMA = EPS*WORK( INDS+R1-1 )
00275       R = R1
00276       DO 110 I = R1, R2 - 1
00277          TMP = WORK( INDS+I ) + WORK( INDP+I )
00278          IF( TMP.EQ.ZERO )
00279      $      TMP = EPS*WORK( INDS+I )
00280          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
00281             MINGMA = TMP
00282             R = I + 1
00283          END IF
00284  110  CONTINUE
00285 *
00286 *     Compute the FP vector: solve N^T v = e_r
00287 *
00288       ISUPPZ( 1 ) = B1
00289       ISUPPZ( 2 ) = BN
00290       Z( R ) = ONE
00291       ZTZ = ONE
00292 *
00293 *     Compute the FP vector upwards from R
00294 *
00295       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00296          DO 210 I = R-1, B1, -1
00297             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00298             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00299      $           THEN
00300                Z( I ) = ZERO
00301                ISUPPZ( 1 ) = I + 1
00302                GOTO 220
00303             ENDIF
00304             ZTZ = ZTZ + Z( I )*Z( I )
00305  210     CONTINUE
00306  220     CONTINUE
00307       ELSE
00308 *        Run slower loop if NaN occurred.
00309          DO 230 I = R - 1, B1, -1
00310             IF( Z( I+1 ).EQ.ZERO ) THEN
00311                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
00312             ELSE
00313                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00314             END IF
00315             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00316      $           THEN
00317                Z( I ) = ZERO
00318                ISUPPZ( 1 ) = I + 1
00319                GO TO 240
00320             END IF
00321             ZTZ = ZTZ + Z( I )*Z( I )
00322  230     CONTINUE
00323  240     CONTINUE
00324       ENDIF
00325 
00326 *     Compute the FP vector downwards from R in blocks of size BLKSIZ
00327       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00328          DO 250 I = R, BN-1
00329             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00330             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00331      $         THEN
00332                Z( I+1 ) = ZERO
00333                ISUPPZ( 2 ) = I
00334                GO TO 260
00335             END IF
00336             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
00337  250     CONTINUE
00338  260     CONTINUE
00339       ELSE
00340 *        Run slower loop if NaN occurred.
00341          DO 270 I = R, BN - 1
00342             IF( Z( I ).EQ.ZERO ) THEN
00343                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
00344             ELSE
00345                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00346             END IF
00347             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00348      $           THEN
00349                Z( I+1 ) = ZERO
00350                ISUPPZ( 2 ) = I
00351                GO TO 280
00352             END IF
00353             ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
00354  270     CONTINUE
00355  280     CONTINUE
00356       END IF
00357 *
00358 *     Compute quantities for convergence test
00359 *
00360       TMP = ONE / ZTZ
00361       NRMINV = SQRT( TMP )
00362       RESID = ABS( MINGMA )*NRMINV
00363       RQCORR = MINGMA*TMP
00364 *
00365 *
00366       RETURN
00367 *
00368 *     End of DLAR1V
00369 *
00370       END