LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlar1v.f
Go to the documentation of this file.
00001 *> \brief \b ZLAR1V
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLAR1V + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlar1v.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlar1v.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlar1v.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
00022 *                  PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
00023 *                  R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       LOGICAL            WANTNC
00027 *       INTEGER   B1, BN, N, NEGCNT, R
00028 *       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
00029 *      $                   RQCORR, ZTZ
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            ISUPPZ( * )
00033 *       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
00034 *      $                  WORK( * )
00035 *       COMPLEX*16       Z( * )
00036 *       ..
00037 *  
00038 *
00039 *> \par Purpose:
00040 *  =============
00041 *>
00042 *> \verbatim
00043 *>
00044 *> ZLAR1V computes the (scaled) r-th column of the inverse of
00045 *> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
00046 *> L D L**T - sigma I. When sigma is close to an eigenvalue, the
00047 *> computed vector is an accurate eigenvector. Usually, r corresponds
00048 *> to the index where the eigenvector is largest in magnitude.
00049 *> The following steps accomplish this computation :
00050 *> (a) Stationary qd transform,  L D L**T - sigma I = L(+) D(+) L(+)**T,
00051 *> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
00052 *> (c) Computation of the diagonal elements of the inverse of
00053 *>     L D L**T - sigma I by combining the above transforms, and choosing
00054 *>     r as the index where the diagonal of the inverse is (one of the)
00055 *>     largest in magnitude.
00056 *> (d) Computation of the (scaled) r-th column of the inverse using the
00057 *>     twisted factorization obtained by combining the top part of the
00058 *>     the stationary and the bottom part of the progressive transform.
00059 *> \endverbatim
00060 *
00061 *  Arguments:
00062 *  ==========
00063 *
00064 *> \param[in] N
00065 *> \verbatim
00066 *>          N is INTEGER
00067 *>           The order of the matrix L D L**T.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] B1
00071 *> \verbatim
00072 *>          B1 is INTEGER
00073 *>           First index of the submatrix of L D L**T.
00074 *> \endverbatim
00075 *>
00076 *> \param[in] BN
00077 *> \verbatim
00078 *>          BN is INTEGER
00079 *>           Last index of the submatrix of L D L**T.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] LAMBDA
00083 *> \verbatim
00084 *>          LAMBDA is DOUBLE PRECISION
00085 *>           The shift. In order to compute an accurate eigenvector,
00086 *>           LAMBDA should be a good approximation to an eigenvalue
00087 *>           of L D L**T.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] L
00091 *> \verbatim
00092 *>          L is DOUBLE PRECISION array, dimension (N-1)
00093 *>           The (n-1) subdiagonal elements of the unit bidiagonal matrix
00094 *>           L, in elements 1 to N-1.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] D
00098 *> \verbatim
00099 *>          D is DOUBLE PRECISION array, dimension (N)
00100 *>           The n diagonal elements of the diagonal matrix D.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] LD
00104 *> \verbatim
00105 *>          LD is DOUBLE PRECISION array, dimension (N-1)
00106 *>           The n-1 elements L(i)*D(i).
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LLD
00110 *> \verbatim
00111 *>          LLD is DOUBLE PRECISION array, dimension (N-1)
00112 *>           The n-1 elements L(i)*L(i)*D(i).
00113 *> \endverbatim
00114 *>
00115 *> \param[in] PIVMIN
00116 *> \verbatim
00117 *>          PIVMIN is DOUBLE PRECISION
00118 *>           The minimum pivot in the Sturm sequence.
00119 *> \endverbatim
00120 *>
00121 *> \param[in] GAPTOL
00122 *> \verbatim
00123 *>          GAPTOL is DOUBLE PRECISION
00124 *>           Tolerance that indicates when eigenvector entries are negligible
00125 *>           w.r.t. their contribution to the residual.
00126 *> \endverbatim
00127 *>
00128 *> \param[in,out] Z
00129 *> \verbatim
00130 *>          Z is COMPLEX*16 array, dimension (N)
00131 *>           On input, all entries of Z must be set to 0.
00132 *>           On output, Z contains the (scaled) r-th column of the
00133 *>           inverse. The scaling is such that Z(R) equals 1.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] WANTNC
00137 *> \verbatim
00138 *>          WANTNC is LOGICAL
00139 *>           Specifies whether NEGCNT has to be computed.
00140 *> \endverbatim
00141 *>
00142 *> \param[out] NEGCNT
00143 *> \verbatim
00144 *>          NEGCNT is INTEGER
00145 *>           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
00146 *>           in the  matrix factorization L D L**T, and NEGCNT = -1 otherwise.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] ZTZ
00150 *> \verbatim
00151 *>          ZTZ is DOUBLE PRECISION
00152 *>           The square of the 2-norm of Z.
00153 *> \endverbatim
00154 *>
00155 *> \param[out] MINGMA
00156 *> \verbatim
00157 *>          MINGMA is DOUBLE PRECISION
00158 *>           The reciprocal of the largest (in magnitude) diagonal
00159 *>           element of the inverse of L D L**T - sigma I.
00160 *> \endverbatim
00161 *>
00162 *> \param[in,out] R
00163 *> \verbatim
00164 *>          R is INTEGER
00165 *>           The twist index for the twisted factorization used to
00166 *>           compute Z.
00167 *>           On input, 0 <= R <= N. If R is input as 0, R is set to
00168 *>           the index where (L D L**T - sigma I)^{-1} is largest
00169 *>           in magnitude. If 1 <= R <= N, R is unchanged.
00170 *>           On output, R contains the twist index used to compute Z.
00171 *>           Ideally, R designates the position of the maximum entry in the
00172 *>           eigenvector.
00173 *> \endverbatim
00174 *>
00175 *> \param[out] ISUPPZ
00176 *> \verbatim
00177 *>          ISUPPZ is INTEGER array, dimension (2)
00178 *>           The support of the vector in Z, i.e., the vector Z is
00179 *>           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
00180 *> \endverbatim
00181 *>
00182 *> \param[out] NRMINV
00183 *> \verbatim
00184 *>          NRMINV is DOUBLE PRECISION
00185 *>           NRMINV = 1/SQRT( ZTZ )
00186 *> \endverbatim
00187 *>
00188 *> \param[out] RESID
00189 *> \verbatim
00190 *>          RESID is DOUBLE PRECISION
00191 *>           The residual of the FP vector.
00192 *>           RESID = ABS( MINGMA )/SQRT( ZTZ )
00193 *> \endverbatim
00194 *>
00195 *> \param[out] RQCORR
00196 *> \verbatim
00197 *>          RQCORR is DOUBLE PRECISION
00198 *>           The Rayleigh Quotient correction to LAMBDA.
00199 *>           RQCORR = MINGMA*TMP
00200 *> \endverbatim
00201 *>
00202 *> \param[out] WORK
00203 *> \verbatim
00204 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
00205 *> \endverbatim
00206 *
00207 *  Authors:
00208 *  ========
00209 *
00210 *> \author Univ. of Tennessee 
00211 *> \author Univ. of California Berkeley 
00212 *> \author Univ. of Colorado Denver 
00213 *> \author NAG Ltd. 
00214 *
00215 *> \date November 2011
00216 *
00217 *> \ingroup complex16OTHERauxiliary
00218 *
00219 *> \par Contributors:
00220 *  ==================
00221 *>
00222 *> Beresford Parlett, University of California, Berkeley, USA \n
00223 *> Jim Demmel, University of California, Berkeley, USA \n
00224 *> Inderjit Dhillon, University of Texas, Austin, USA \n
00225 *> Osni Marques, LBNL/NERSC, USA \n
00226 *> Christof Voemel, University of California, Berkeley, USA
00227 *
00228 *  =====================================================================
00229       SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
00230      $           PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
00231      $           R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
00232 *
00233 *  -- LAPACK auxiliary routine (version 3.4.0) --
00234 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00235 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00236 *     November 2011
00237 *
00238 *     .. Scalar Arguments ..
00239       LOGICAL            WANTNC
00240       INTEGER   B1, BN, N, NEGCNT, R
00241       DOUBLE PRECISION   GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
00242      $                   RQCORR, ZTZ
00243 *     ..
00244 *     .. Array Arguments ..
00245       INTEGER            ISUPPZ( * )
00246       DOUBLE PRECISION   D( * ), L( * ), LD( * ), LLD( * ),
00247      $                  WORK( * )
00248       COMPLEX*16       Z( * )
00249 *     ..
00250 *
00251 *  =====================================================================
00252 *
00253 *     .. Parameters ..
00254       DOUBLE PRECISION   ZERO, ONE
00255       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00256       COMPLEX*16         CONE
00257       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
00258 
00259 *     ..
00260 *     .. Local Scalars ..
00261       LOGICAL            SAWNAN1, SAWNAN2
00262       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
00263      $                   R2
00264       DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
00265 *     ..
00266 *     .. External Functions ..
00267       LOGICAL DISNAN
00268       DOUBLE PRECISION   DLAMCH
00269       EXTERNAL           DISNAN, DLAMCH
00270 *     ..
00271 *     .. Intrinsic Functions ..
00272       INTRINSIC          ABS, DBLE
00273 *     ..
00274 *     .. Executable Statements ..
00275 *
00276       EPS = DLAMCH( 'Precision' )
00277 
00278 
00279       IF( R.EQ.0 ) THEN
00280          R1 = B1
00281          R2 = BN
00282       ELSE
00283          R1 = R
00284          R2 = R
00285       END IF
00286 
00287 *     Storage for LPLUS
00288       INDLPL = 0
00289 *     Storage for UMINUS
00290       INDUMN = N
00291       INDS = 2*N + 1
00292       INDP = 3*N + 1
00293 
00294       IF( B1.EQ.1 ) THEN
00295          WORK( INDS ) = ZERO
00296       ELSE
00297          WORK( INDS+B1-1 ) = LLD( B1-1 )
00298       END IF
00299 
00300 *
00301 *     Compute the stationary transform (using the differential form)
00302 *     until the index R2.
00303 *
00304       SAWNAN1 = .FALSE.
00305       NEG1 = 0
00306       S = WORK( INDS+B1-1 ) - LAMBDA
00307       DO 50 I = B1, R1 - 1
00308          DPLUS = D( I ) + S
00309          WORK( INDLPL+I ) = LD( I ) / DPLUS
00310          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00311          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00312          S = WORK( INDS+I ) - LAMBDA
00313  50   CONTINUE
00314       SAWNAN1 = DISNAN( S )
00315       IF( SAWNAN1 ) GOTO 60
00316       DO 51 I = R1, R2 - 1
00317          DPLUS = D( I ) + S
00318          WORK( INDLPL+I ) = LD( I ) / DPLUS
00319          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00320          S = WORK( INDS+I ) - LAMBDA
00321  51   CONTINUE
00322       SAWNAN1 = DISNAN( S )
00323 *
00324  60   CONTINUE
00325       IF( SAWNAN1 ) THEN
00326 *        Runs a slower version of the above loop if a NaN is detected
00327          NEG1 = 0
00328          S = WORK( INDS+B1-1 ) - LAMBDA
00329          DO 70 I = B1, R1 - 1
00330             DPLUS = D( I ) + S
00331             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00332             WORK( INDLPL+I ) = LD( I ) / DPLUS
00333             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00334             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00335             IF( WORK( INDLPL+I ).EQ.ZERO )
00336      $                      WORK( INDS+I ) = LLD( I )
00337             S = WORK( INDS+I ) - LAMBDA
00338  70      CONTINUE
00339          DO 71 I = R1, R2 - 1
00340             DPLUS = D( I ) + S
00341             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00342             WORK( INDLPL+I ) = LD( I ) / DPLUS
00343             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00344             IF( WORK( INDLPL+I ).EQ.ZERO )
00345      $                      WORK( INDS+I ) = LLD( I )
00346             S = WORK( INDS+I ) - LAMBDA
00347  71      CONTINUE
00348       END IF
00349 *
00350 *     Compute the progressive transform (using the differential form)
00351 *     until the index R1
00352 *
00353       SAWNAN2 = .FALSE.
00354       NEG2 = 0
00355       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
00356       DO 80 I = BN - 1, R1, -1
00357          DMINUS = LLD( I ) + WORK( INDP+I )
00358          TMP = D( I ) / DMINUS
00359          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00360          WORK( INDUMN+I ) = L( I )*TMP
00361          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00362  80   CONTINUE
00363       TMP = WORK( INDP+R1-1 )
00364       SAWNAN2 = DISNAN( TMP )
00365 
00366       IF( SAWNAN2 ) THEN
00367 *        Runs a slower version of the above loop if a NaN is detected
00368          NEG2 = 0
00369          DO 100 I = BN-1, R1, -1
00370             DMINUS = LLD( I ) + WORK( INDP+I )
00371             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
00372             TMP = D( I ) / DMINUS
00373             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00374             WORK( INDUMN+I ) = L( I )*TMP
00375             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00376             IF( TMP.EQ.ZERO )
00377      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
00378  100     CONTINUE
00379       END IF
00380 *
00381 *     Find the index (from R1 to R2) of the largest (in magnitude)
00382 *     diagonal element of the inverse
00383 *
00384       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
00385       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
00386       IF( WANTNC ) THEN
00387          NEGCNT = NEG1 + NEG2
00388       ELSE
00389          NEGCNT = -1
00390       ENDIF
00391       IF( ABS(MINGMA).EQ.ZERO )
00392      $   MINGMA = EPS*WORK( INDS+R1-1 )
00393       R = R1
00394       DO 110 I = R1, R2 - 1
00395          TMP = WORK( INDS+I ) + WORK( INDP+I )
00396          IF( TMP.EQ.ZERO )
00397      $      TMP = EPS*WORK( INDS+I )
00398          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
00399             MINGMA = TMP
00400             R = I + 1
00401          END IF
00402  110  CONTINUE
00403 *
00404 *     Compute the FP vector: solve N^T v = e_r
00405 *
00406       ISUPPZ( 1 ) = B1
00407       ISUPPZ( 2 ) = BN
00408       Z( R ) = CONE
00409       ZTZ = ONE
00410 *
00411 *     Compute the FP vector upwards from R
00412 *
00413       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00414          DO 210 I = R-1, B1, -1
00415             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00416             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00417      $           THEN
00418                Z( I ) = ZERO
00419                ISUPPZ( 1 ) = I + 1
00420                GOTO 220
00421             ENDIF
00422             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
00423  210     CONTINUE
00424  220     CONTINUE
00425       ELSE
00426 *        Run slower loop if NaN occurred.
00427          DO 230 I = R - 1, B1, -1
00428             IF( Z( I+1 ).EQ.ZERO ) THEN
00429                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
00430             ELSE
00431                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00432             END IF
00433             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00434      $           THEN
00435                Z( I ) = ZERO
00436                ISUPPZ( 1 ) = I + 1
00437                GO TO 240
00438             END IF
00439             ZTZ = ZTZ + DBLE( Z( I )*Z( I ) )
00440  230     CONTINUE
00441  240     CONTINUE
00442       ENDIF
00443 
00444 *     Compute the FP vector downwards from R in blocks of size BLKSIZ
00445       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00446          DO 250 I = R, BN-1
00447             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00448             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00449      $         THEN
00450                Z( I+1 ) = ZERO
00451                ISUPPZ( 2 ) = I
00452                GO TO 260
00453             END IF
00454             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
00455  250     CONTINUE
00456  260     CONTINUE
00457       ELSE
00458 *        Run slower loop if NaN occurred.
00459          DO 270 I = R, BN - 1
00460             IF( Z( I ).EQ.ZERO ) THEN
00461                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
00462             ELSE
00463                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00464             END IF
00465             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00466      $           THEN
00467                Z( I+1 ) = ZERO
00468                ISUPPZ( 2 ) = I
00469                GO TO 280
00470             END IF
00471             ZTZ = ZTZ + DBLE( Z( I+1 )*Z( I+1 ) )
00472  270     CONTINUE
00473  280     CONTINUE
00474       END IF
00475 *
00476 *     Compute quantities for convergence test
00477 *
00478       TMP = ONE / ZTZ
00479       NRMINV = SQRT( TMP )
00480       RESID = ABS( MINGMA )*NRMINV
00481       RQCORR = MINGMA*TMP
00482 *
00483 *
00484       RETURN
00485 *
00486 *     End of ZLAR1V
00487 *
00488       END
 All Files Functions