LAPACK: SRC/clar1v.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAR1V( 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       REAL               GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
00014      $                   RQCORR, ZTZ
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            ISUPPZ( * )
00018       REAL               D( * ), L( * ), LD( * ), LLD( * ),
00019      $                  WORK( * )
00020       COMPLEX          Z( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  CLAR1V 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) REAL            
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) REAL             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) REAL             array, dimension (N)
00064 *           The n diagonal elements of the diagonal matrix D.
00065 *
00066 *  LD       (input) REAL             array, dimension (N-1)
00067 *           The n-1 elements L(i)*D(i).
00068 *
00069 *  LLD      (input) REAL             array, dimension (N-1)
00070 *           The n-1 elements L(i)*L(i)*D(i).
00071 *
00072 *  PIVMIN   (input) REAL            
00073 *           The minimum pivot in the Sturm sequence.
00074 *
00075 *  GAPTOL   (input) REAL            
00076 *           Tolerance that indicates when eigenvector entries are negligible
00077 *           w.r.t. their contribution to the residual.
00078 *
00079 *  Z        (input/output) COMPLEX          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) REAL            
00092 *           The square of the 2-norm of Z.
00093 *
00094 *  MINGMA   (output) REAL            
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) REAL            
00113 *           NRMINV = 1/SQRT( ZTZ )
00114 *
00115 *  RESID    (output) REAL            
00116 *           The residual of the FP vector.
00117 *           RESID = ABS( MINGMA )/SQRT( ZTZ )
00118 *
00119 *  RQCORR   (output) REAL            
00120 *           The Rayleigh Quotient correction to LAMBDA.
00121 *           RQCORR = MINGMA*TMP
00122 *
00123 *  WORK     (workspace) REAL             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       REAL               ZERO, ONE
00139       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00140       COMPLEX            CONE
00141       PARAMETER          ( CONE = ( 1.0E0, 0.0E0 ) )
00142 
00143 *     ..
00144 *     .. Local Scalars ..
00145       LOGICAL            SAWNAN1, SAWNAN2
00146       INTEGER            I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
00147      $                   R2
00148       REAL               DMINUS, DPLUS, EPS, S, TMP
00149 *     ..
00150 *     .. External Functions ..
00151       LOGICAL SISNAN
00152       REAL               SLAMCH
00153       EXTERNAL           SISNAN, SLAMCH
00154 *     ..
00155 *     .. Intrinsic Functions ..
00156       INTRINSIC          ABS, REAL
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160       EPS = SLAMCH( 'Precision' )
00161 
00162 
00163       IF( R.EQ.0 ) THEN
00164          R1 = B1
00165          R2 = BN
00166       ELSE
00167          R1 = R
00168          R2 = R
00169       END IF
00170 
00171 *     Storage for LPLUS
00172       INDLPL = 0
00173 *     Storage for UMINUS
00174       INDUMN = N
00175       INDS = 2*N + 1
00176       INDP = 3*N + 1
00177 
00178       IF( B1.EQ.1 ) THEN
00179          WORK( INDS ) = ZERO
00180       ELSE
00181          WORK( INDS+B1-1 ) = LLD( B1-1 )
00182       END IF
00183 
00184 *
00185 *     Compute the stationary transform (using the differential form)
00186 *     until the index R2.
00187 *
00188       SAWNAN1 = .FALSE.
00189       NEG1 = 0
00190       S = WORK( INDS+B1-1 ) - LAMBDA
00191       DO 50 I = B1, R1 - 1
00192          DPLUS = D( I ) + S
00193          WORK( INDLPL+I ) = LD( I ) / DPLUS
00194          IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00195          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00196          S = WORK( INDS+I ) - LAMBDA
00197  50   CONTINUE
00198       SAWNAN1 = SISNAN( S )
00199       IF( SAWNAN1 ) GOTO 60
00200       DO 51 I = R1, R2 - 1
00201          DPLUS = D( I ) + S
00202          WORK( INDLPL+I ) = LD( I ) / DPLUS
00203          WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00204          S = WORK( INDS+I ) - LAMBDA
00205  51   CONTINUE
00206       SAWNAN1 = SISNAN( S )
00207 *
00208  60   CONTINUE
00209       IF( SAWNAN1 ) THEN
00210 *        Runs a slower version of the above loop if a NaN is detected
00211          NEG1 = 0
00212          S = WORK( INDS+B1-1 ) - LAMBDA
00213          DO 70 I = B1, R1 - 1
00214             DPLUS = D( I ) + S
00215             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00216             WORK( INDLPL+I ) = LD( I ) / DPLUS
00217             IF(DPLUS.LT.ZERO) NEG1 = NEG1 + 1
00218             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00219             IF( WORK( INDLPL+I ).EQ.ZERO )
00220      $                      WORK( INDS+I ) = LLD( I )
00221             S = WORK( INDS+I ) - LAMBDA
00222  70      CONTINUE
00223          DO 71 I = R1, R2 - 1
00224             DPLUS = D( I ) + S
00225             IF(ABS(DPLUS).LT.PIVMIN) DPLUS = -PIVMIN
00226             WORK( INDLPL+I ) = LD( I ) / DPLUS
00227             WORK( INDS+I ) = S*WORK( INDLPL+I )*L( I )
00228             IF( WORK( INDLPL+I ).EQ.ZERO )
00229      $                      WORK( INDS+I ) = LLD( I )
00230             S = WORK( INDS+I ) - LAMBDA
00231  71      CONTINUE
00232       END IF
00233 *
00234 *     Compute the progressive transform (using the differential form)
00235 *     until the index R1
00236 *
00237       SAWNAN2 = .FALSE.
00238       NEG2 = 0
00239       WORK( INDP+BN-1 ) = D( BN ) - LAMBDA
00240       DO 80 I = BN - 1, R1, -1
00241          DMINUS = LLD( I ) + WORK( INDP+I )
00242          TMP = D( I ) / DMINUS
00243          IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00244          WORK( INDUMN+I ) = L( I )*TMP
00245          WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00246  80   CONTINUE
00247       TMP = WORK( INDP+R1-1 )
00248       SAWNAN2 = SISNAN( TMP )
00249 
00250       IF( SAWNAN2 ) THEN
00251 *        Runs a slower version of the above loop if a NaN is detected
00252          NEG2 = 0
00253          DO 100 I = BN-1, R1, -1
00254             DMINUS = LLD( I ) + WORK( INDP+I )
00255             IF(ABS(DMINUS).LT.PIVMIN) DMINUS = -PIVMIN
00256             TMP = D( I ) / DMINUS
00257             IF(DMINUS.LT.ZERO) NEG2 = NEG2 + 1
00258             WORK( INDUMN+I ) = L( I )*TMP
00259             WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - LAMBDA
00260             IF( TMP.EQ.ZERO )
00261      $          WORK( INDP+I-1 ) = D( I ) - LAMBDA
00262  100     CONTINUE
00263       END IF
00264 *
00265 *     Find the index (from R1 to R2) of the largest (in magnitude)
00266 *     diagonal element of the inverse
00267 *
00268       MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
00269       IF( MINGMA.LT.ZERO ) NEG1 = NEG1 + 1
00270       IF( WANTNC ) THEN
00271          NEGCNT = NEG1 + NEG2
00272       ELSE
00273          NEGCNT = -1
00274       ENDIF
00275       IF( ABS(MINGMA).EQ.ZERO )
00276      $   MINGMA = EPS*WORK( INDS+R1-1 )
00277       R = R1
00278       DO 110 I = R1, R2 - 1
00279          TMP = WORK( INDS+I ) + WORK( INDP+I )
00280          IF( TMP.EQ.ZERO )
00281      $      TMP = EPS*WORK( INDS+I )
00282          IF( ABS( TMP ).LE.ABS( MINGMA ) ) THEN
00283             MINGMA = TMP
00284             R = I + 1
00285          END IF
00286  110  CONTINUE
00287 *
00288 *     Compute the FP vector: solve N^T v = e_r
00289 *
00290       ISUPPZ( 1 ) = B1
00291       ISUPPZ( 2 ) = BN
00292       Z( R ) = CONE
00293       ZTZ = ONE
00294 *
00295 *     Compute the FP vector upwards from R
00296 *
00297       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00298          DO 210 I = R-1, B1, -1
00299             Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00300             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00301      $           THEN
00302                Z( I ) = ZERO
00303                ISUPPZ( 1 ) = I + 1
00304                GOTO 220
00305             ENDIF
00306             ZTZ = ZTZ + REAL( Z( I )*Z( I ) )
00307  210     CONTINUE
00308  220     CONTINUE
00309       ELSE
00310 *        Run slower loop if NaN occurred.
00311          DO 230 I = R - 1, B1, -1
00312             IF( Z( I+1 ).EQ.ZERO ) THEN
00313                Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
00314             ELSE
00315                Z( I ) = -( WORK( INDLPL+I )*Z( I+1 ) )
00316             END IF
00317             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00318      $           THEN
00319                Z( I ) = ZERO
00320                ISUPPZ( 1 ) = I + 1
00321                GO TO 240
00322             END IF
00323             ZTZ = ZTZ + REAL( Z( I )*Z( I ) )
00324  230     CONTINUE
00325  240     CONTINUE
00326       ENDIF
00327 
00328 *     Compute the FP vector downwards from R in blocks of size BLKSIZ
00329       IF( .NOT.SAWNAN1 .AND. .NOT.SAWNAN2 ) THEN
00330          DO 250 I = R, BN-1
00331             Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00332             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00333      $         THEN
00334                Z( I+1 ) = ZERO
00335                ISUPPZ( 2 ) = I
00336                GO TO 260
00337             END IF
00338             ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) )
00339  250     CONTINUE
00340  260     CONTINUE
00341       ELSE
00342 *        Run slower loop if NaN occurred.
00343          DO 270 I = R, BN - 1
00344             IF( Z( I ).EQ.ZERO ) THEN
00345                Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
00346             ELSE
00347                Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
00348             END IF
00349             IF( (ABS(Z(I))+ABS(Z(I+1)))* ABS(LD(I)).LT.GAPTOL )
00350      $           THEN
00351                Z( I+1 ) = ZERO
00352                ISUPPZ( 2 ) = I
00353                GO TO 280
00354             END IF
00355             ZTZ = ZTZ + REAL( Z( I+1 )*Z( I+1 ) )
00356  270     CONTINUE
00357  280     CONTINUE
00358       END IF
00359 *
00360 *     Compute quantities for convergence test
00361 *
00362       TMP = ONE / ZTZ
00363       NRMINV = SQRT( TMP )
00364       RESID = ABS( MINGMA )*NRMINV
00365       RQCORR = MINGMA*TMP
00366 *
00367 *
00368       RETURN
00369 *
00370 *     End of CLAR1V
00371 *
00372       END