LAPACK 3.3.0

zhsein.f

Go to the documentation of this file.
00001       SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
00002      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
00003      $                   IFAILR, INFO )
00004 *
00005 *  -- LAPACK 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       CHARACTER          EIGSRC, INITV, SIDE
00012       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            SELECT( * )
00016       INTEGER            IFAILL( * ), IFAILR( * )
00017       DOUBLE PRECISION   RWORK( * )
00018       COMPLEX*16         H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00019      $                   W( * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  ZHSEIN uses inverse iteration to find specified right and/or left
00026 *  eigenvectors of a complex upper Hessenberg matrix H.
00027 *
00028 *  The right eigenvector x and the left eigenvector y of the matrix H
00029 *  corresponding to an eigenvalue w are defined by:
00030 *
00031 *               H * x = w * x,     y**h * H = w * y**h
00032 *
00033 *  where y**h denotes the conjugate transpose of the vector y.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  SIDE    (input) CHARACTER*1
00039 *          = 'R': compute right eigenvectors only;
00040 *          = 'L': compute left eigenvectors only;
00041 *          = 'B': compute both right and left eigenvectors.
00042 *
00043 *  EIGSRC  (input) CHARACTER*1
00044 *          Specifies the source of eigenvalues supplied in W:
00045 *          = 'Q': the eigenvalues were found using ZHSEQR; thus, if
00046 *                 H has zero subdiagonal elements, and so is
00047 *                 block-triangular, then the j-th eigenvalue can be
00048 *                 assumed to be an eigenvalue of the block containing
00049 *                 the j-th row/column.  This property allows ZHSEIN to
00050 *                 perform inverse iteration on just one diagonal block.
00051 *          = 'N': no assumptions are made on the correspondence
00052 *                 between eigenvalues and diagonal blocks.  In this
00053 *                 case, ZHSEIN must always perform inverse iteration
00054 *                 using the whole matrix H.
00055 *
00056 *  INITV   (input) CHARACTER*1
00057 *          = 'N': no initial vectors are supplied;
00058 *          = 'U': user-supplied initial vectors are stored in the arrays
00059 *                 VL and/or VR.
00060 *
00061 *  SELECT  (input) LOGICAL array, dimension (N)
00062 *          Specifies the eigenvectors to be computed. To select the
00063 *          eigenvector corresponding to the eigenvalue W(j),
00064 *          SELECT(j) must be set to .TRUE..
00065 *
00066 *  N       (input) INTEGER
00067 *          The order of the matrix H.  N >= 0.
00068 *
00069 *  H       (input) COMPLEX*16 array, dimension (LDH,N)
00070 *          The upper Hessenberg matrix H.
00071 *
00072 *  LDH     (input) INTEGER
00073 *          The leading dimension of the array H.  LDH >= max(1,N).
00074 *
00075 *  W       (input/output) COMPLEX*16 array, dimension (N)
00076 *          On entry, the eigenvalues of H.
00077 *          On exit, the real parts of W may have been altered since
00078 *          close eigenvalues are perturbed slightly in searching for
00079 *          independent eigenvectors.
00080 *
00081 *  VL      (input/output) COMPLEX*16 array, dimension (LDVL,MM)
00082 *          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
00083 *          contain starting vectors for the inverse iteration for the
00084 *          left eigenvectors; the starting vector for each eigenvector
00085 *          must be in the same column in which the eigenvector will be
00086 *          stored.
00087 *          On exit, if SIDE = 'L' or 'B', the left eigenvectors
00088 *          specified by SELECT will be stored consecutively in the
00089 *          columns of VL, in the same order as their eigenvalues.
00090 *          If SIDE = 'R', VL is not referenced.
00091 *
00092 *  LDVL    (input) INTEGER
00093 *          The leading dimension of the array VL.
00094 *          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
00095 *
00096 *  VR      (input/output) COMPLEX*16 array, dimension (LDVR,MM)
00097 *          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
00098 *          contain starting vectors for the inverse iteration for the
00099 *          right eigenvectors; the starting vector for each eigenvector
00100 *          must be in the same column in which the eigenvector will be
00101 *          stored.
00102 *          On exit, if SIDE = 'R' or 'B', the right eigenvectors
00103 *          specified by SELECT will be stored consecutively in the
00104 *          columns of VR, in the same order as their eigenvalues.
00105 *          If SIDE = 'L', VR is not referenced.
00106 *
00107 *  LDVR    (input) INTEGER
00108 *          The leading dimension of the array VR.
00109 *          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
00110 *
00111 *  MM      (input) INTEGER
00112 *          The number of columns in the arrays VL and/or VR. MM >= M.
00113 *
00114 *  M       (output) INTEGER
00115 *          The number of columns in the arrays VL and/or VR required to
00116 *          store the eigenvectors (= the number of .TRUE. elements in
00117 *          SELECT).
00118 *
00119 *  WORK    (workspace) COMPLEX*16 array, dimension (N*N)
00120 *
00121 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00122 *
00123 *  IFAILL  (output) INTEGER array, dimension (MM)
00124 *          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
00125 *          eigenvector in the i-th column of VL (corresponding to the
00126 *          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
00127 *          eigenvector converged satisfactorily.
00128 *          If SIDE = 'R', IFAILL is not referenced.
00129 *
00130 *  IFAILR  (output) INTEGER array, dimension (MM)
00131 *          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
00132 *          eigenvector in the i-th column of VR (corresponding to the
00133 *          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
00134 *          eigenvector converged satisfactorily.
00135 *          If SIDE = 'L', IFAILR is not referenced.
00136 *
00137 *  INFO    (output) INTEGER
00138 *          = 0:  successful exit
00139 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00140 *          > 0:  if INFO = i, i is the number of eigenvectors which
00141 *                failed to converge; see IFAILL and IFAILR for further
00142 *                details.
00143 *
00144 *  Further Details
00145 *  ===============
00146 *
00147 *  Each eigenvector is normalized so that the element of largest
00148 *  magnitude has magnitude 1; here the magnitude of a complex number
00149 *  (x,y) is taken to be |x|+|y|.
00150 *
00151 *  =====================================================================
00152 *
00153 *     .. Parameters ..
00154       COMPLEX*16         ZERO
00155       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
00156       DOUBLE PRECISION   RZERO
00157       PARAMETER          ( RZERO = 0.0D+0 )
00158 *     ..
00159 *     .. Local Scalars ..
00160       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
00161       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
00162       DOUBLE PRECISION   EPS3, HNORM, SMLNUM, ULP, UNFL
00163       COMPLEX*16         CDUM, WK
00164 *     ..
00165 *     .. External Functions ..
00166       LOGICAL            LSAME
00167       DOUBLE PRECISION   DLAMCH, ZLANHS
00168       EXTERNAL           LSAME, DLAMCH, ZLANHS
00169 *     ..
00170 *     .. External Subroutines ..
00171       EXTERNAL           XERBLA, ZLAEIN
00172 *     ..
00173 *     .. Intrinsic Functions ..
00174       INTRINSIC          ABS, DBLE, DIMAG, MAX
00175 *     ..
00176 *     .. Statement Functions ..
00177       DOUBLE PRECISION   CABS1
00178 *     ..
00179 *     .. Statement Function definitions ..
00180       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
00181 *     ..
00182 *     .. Executable Statements ..
00183 *
00184 *     Decode and test the input parameters.
00185 *
00186       BOTHV = LSAME( SIDE, 'B' )
00187       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00188       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00189 *
00190       FROMQR = LSAME( EIGSRC, 'Q' )
00191 *
00192       NOINIT = LSAME( INITV, 'N' )
00193 *
00194 *     Set M to the number of columns required to store the selected
00195 *     eigenvectors.
00196 *
00197       M = 0
00198       DO 10 K = 1, N
00199          IF( SELECT( K ) )
00200      $      M = M + 1
00201    10 CONTINUE
00202 *
00203       INFO = 0
00204       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00205          INFO = -1
00206       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
00207          INFO = -2
00208       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
00209          INFO = -3
00210       ELSE IF( N.LT.0 ) THEN
00211          INFO = -5
00212       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
00213          INFO = -7
00214       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00215          INFO = -10
00216       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00217          INFO = -12
00218       ELSE IF( MM.LT.M ) THEN
00219          INFO = -13
00220       END IF
00221       IF( INFO.NE.0 ) THEN
00222          CALL XERBLA( 'ZHSEIN', -INFO )
00223          RETURN
00224       END IF
00225 *
00226 *     Quick return if possible.
00227 *
00228       IF( N.EQ.0 )
00229      $   RETURN
00230 *
00231 *     Set machine-dependent constants.
00232 *
00233       UNFL = DLAMCH( 'Safe minimum' )
00234       ULP = DLAMCH( 'Precision' )
00235       SMLNUM = UNFL*( N / ULP )
00236 *
00237       LDWORK = N
00238 *
00239       KL = 1
00240       KLN = 0
00241       IF( FROMQR ) THEN
00242          KR = 0
00243       ELSE
00244          KR = N
00245       END IF
00246       KS = 1
00247 *
00248       DO 100 K = 1, N
00249          IF( SELECT( K ) ) THEN
00250 *
00251 *           Compute eigenvector(s) corresponding to W(K).
00252 *
00253             IF( FROMQR ) THEN
00254 *
00255 *              If affiliation of eigenvalues is known, check whether
00256 *              the matrix splits.
00257 *
00258 *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
00259 *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
00260 *              KR = N).
00261 *
00262 *              Then inverse iteration can be performed with the
00263 *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
00264 *              the submatrix H(1:KR,1:KR) for a right eigenvector.
00265 *
00266                DO 20 I = K, KL + 1, -1
00267                   IF( H( I, I-1 ).EQ.ZERO )
00268      $               GO TO 30
00269    20          CONTINUE
00270    30          CONTINUE
00271                KL = I
00272                IF( K.GT.KR ) THEN
00273                   DO 40 I = K, N - 1
00274                      IF( H( I+1, I ).EQ.ZERO )
00275      $                  GO TO 50
00276    40             CONTINUE
00277    50             CONTINUE
00278                   KR = I
00279                END IF
00280             END IF
00281 *
00282             IF( KL.NE.KLN ) THEN
00283                KLN = KL
00284 *
00285 *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
00286 *              has not ben computed before.
00287 *
00288                HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
00289                IF( HNORM.GT.RZERO ) THEN
00290                   EPS3 = HNORM*ULP
00291                ELSE
00292                   EPS3 = SMLNUM
00293                END IF
00294             END IF
00295 *
00296 *           Perturb eigenvalue if it is close to any previous
00297 *           selected eigenvalues affiliated to the submatrix
00298 *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
00299 *
00300             WK = W( K )
00301    60       CONTINUE
00302             DO 70 I = K - 1, KL, -1
00303                IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
00304                   WK = WK + EPS3
00305                   GO TO 60
00306                END IF
00307    70       CONTINUE
00308             W( K ) = WK
00309 *
00310             IF( LEFTV ) THEN
00311 *
00312 *              Compute left eigenvector.
00313 *
00314                CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
00315      $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
00316      $                      SMLNUM, IINFO )
00317                IF( IINFO.GT.0 ) THEN
00318                   INFO = INFO + 1
00319                   IFAILL( KS ) = K
00320                ELSE
00321                   IFAILL( KS ) = 0
00322                END IF
00323                DO 80 I = 1, KL - 1
00324                   VL( I, KS ) = ZERO
00325    80          CONTINUE
00326             END IF
00327             IF( RIGHTV ) THEN
00328 *
00329 *              Compute right eigenvector.
00330 *
00331                CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
00332      $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
00333                IF( IINFO.GT.0 ) THEN
00334                   INFO = INFO + 1
00335                   IFAILR( KS ) = K
00336                ELSE
00337                   IFAILR( KS ) = 0
00338                END IF
00339                DO 90 I = KR + 1, N
00340                   VR( I, KS ) = ZERO
00341    90          CONTINUE
00342             END IF
00343             KS = KS + 1
00344          END IF
00345   100 CONTINUE
00346 *
00347       RETURN
00348 *
00349 *     End of ZHSEIN
00350 *
00351       END
 All Files Functions