LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhsein.f
Go to the documentation of this file.
00001 *> \brief \b ZHSEIN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHSEIN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhsein.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhsein.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhsein.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
00022 *                          LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
00023 *                          IFAILR, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          EIGSRC, INITV, SIDE
00027 *       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       LOGICAL            SELECT( * )
00031 *       INTEGER            IFAILL( * ), IFAILR( * )
00032 *       DOUBLE PRECISION   RWORK( * )
00033 *       COMPLEX*16         H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00034 *      $                   W( * ), WORK( * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> ZHSEIN uses inverse iteration to find specified right and/or left
00044 *> eigenvectors of a complex upper Hessenberg matrix H.
00045 *>
00046 *> The right eigenvector x and the left eigenvector y of the matrix H
00047 *> corresponding to an eigenvalue w are defined by:
00048 *>
00049 *>              H * x = w * x,     y**h * H = w * y**h
00050 *>
00051 *> where y**h denotes the conjugate transpose of the vector y.
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] SIDE
00058 *> \verbatim
00059 *>          SIDE is CHARACTER*1
00060 *>          = 'R': compute right eigenvectors only;
00061 *>          = 'L': compute left eigenvectors only;
00062 *>          = 'B': compute both right and left eigenvectors.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] EIGSRC
00066 *> \verbatim
00067 *>          EIGSRC is CHARACTER*1
00068 *>          Specifies the source of eigenvalues supplied in W:
00069 *>          = 'Q': the eigenvalues were found using ZHSEQR; thus, if
00070 *>                 H has zero subdiagonal elements, and so is
00071 *>                 block-triangular, then the j-th eigenvalue can be
00072 *>                 assumed to be an eigenvalue of the block containing
00073 *>                 the j-th row/column.  This property allows ZHSEIN to
00074 *>                 perform inverse iteration on just one diagonal block.
00075 *>          = 'N': no assumptions are made on the correspondence
00076 *>                 between eigenvalues and diagonal blocks.  In this
00077 *>                 case, ZHSEIN must always perform inverse iteration
00078 *>                 using the whole matrix H.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] INITV
00082 *> \verbatim
00083 *>          INITV is CHARACTER*1
00084 *>          = 'N': no initial vectors are supplied;
00085 *>          = 'U': user-supplied initial vectors are stored in the arrays
00086 *>                 VL and/or VR.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] SELECT
00090 *> \verbatim
00091 *>          SELECT is LOGICAL array, dimension (N)
00092 *>          Specifies the eigenvectors to be computed. To select the
00093 *>          eigenvector corresponding to the eigenvalue W(j),
00094 *>          SELECT(j) must be set to .TRUE..
00095 *> \endverbatim
00096 *>
00097 *> \param[in] N
00098 *> \verbatim
00099 *>          N is INTEGER
00100 *>          The order of the matrix H.  N >= 0.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] H
00104 *> \verbatim
00105 *>          H is COMPLEX*16 array, dimension (LDH,N)
00106 *>          The upper Hessenberg matrix H.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDH
00110 *> \verbatim
00111 *>          LDH is INTEGER
00112 *>          The leading dimension of the array H.  LDH >= max(1,N).
00113 *> \endverbatim
00114 *>
00115 *> \param[in,out] W
00116 *> \verbatim
00117 *>          W is COMPLEX*16 array, dimension (N)
00118 *>          On entry, the eigenvalues of H.
00119 *>          On exit, the real parts of W may have been altered since
00120 *>          close eigenvalues are perturbed slightly in searching for
00121 *>          independent eigenvectors.
00122 *> \endverbatim
00123 *>
00124 *> \param[in,out] VL
00125 *> \verbatim
00126 *>          VL is COMPLEX*16 array, dimension (LDVL,MM)
00127 *>          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
00128 *>          contain starting vectors for the inverse iteration for the
00129 *>          left eigenvectors; the starting vector for each eigenvector
00130 *>          must be in the same column in which the eigenvector will be
00131 *>          stored.
00132 *>          On exit, if SIDE = 'L' or 'B', the left eigenvectors
00133 *>          specified by SELECT will be stored consecutively in the
00134 *>          columns of VL, in the same order as their eigenvalues.
00135 *>          If SIDE = 'R', VL is not referenced.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] LDVL
00139 *> \verbatim
00140 *>          LDVL is INTEGER
00141 *>          The leading dimension of the array VL.
00142 *>          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
00143 *> \endverbatim
00144 *>
00145 *> \param[in,out] VR
00146 *> \verbatim
00147 *>          VR is COMPLEX*16 array, dimension (LDVR,MM)
00148 *>          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
00149 *>          contain starting vectors for the inverse iteration for the
00150 *>          right eigenvectors; the starting vector for each eigenvector
00151 *>          must be in the same column in which the eigenvector will be
00152 *>          stored.
00153 *>          On exit, if SIDE = 'R' or 'B', the right eigenvectors
00154 *>          specified by SELECT will be stored consecutively in the
00155 *>          columns of VR, in the same order as their eigenvalues.
00156 *>          If SIDE = 'L', VR is not referenced.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] LDVR
00160 *> \verbatim
00161 *>          LDVR is INTEGER
00162 *>          The leading dimension of the array VR.
00163 *>          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
00164 *> \endverbatim
00165 *>
00166 *> \param[in] MM
00167 *> \verbatim
00168 *>          MM is INTEGER
00169 *>          The number of columns in the arrays VL and/or VR. MM >= M.
00170 *> \endverbatim
00171 *>
00172 *> \param[out] M
00173 *> \verbatim
00174 *>          M is INTEGER
00175 *>          The number of columns in the arrays VL and/or VR required to
00176 *>          store the eigenvectors (= the number of .TRUE. elements in
00177 *>          SELECT).
00178 *> \endverbatim
00179 *>
00180 *> \param[out] WORK
00181 *> \verbatim
00182 *>          WORK is COMPLEX*16 array, dimension (N*N)
00183 *> \endverbatim
00184 *>
00185 *> \param[out] RWORK
00186 *> \verbatim
00187 *>          RWORK is DOUBLE PRECISION array, dimension (N)
00188 *> \endverbatim
00189 *>
00190 *> \param[out] IFAILL
00191 *> \verbatim
00192 *>          IFAILL is INTEGER array, dimension (MM)
00193 *>          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
00194 *>          eigenvector in the i-th column of VL (corresponding to the
00195 *>          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
00196 *>          eigenvector converged satisfactorily.
00197 *>          If SIDE = 'R', IFAILL is not referenced.
00198 *> \endverbatim
00199 *>
00200 *> \param[out] IFAILR
00201 *> \verbatim
00202 *>          IFAILR is INTEGER array, dimension (MM)
00203 *>          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
00204 *>          eigenvector in the i-th column of VR (corresponding to the
00205 *>          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
00206 *>          eigenvector converged satisfactorily.
00207 *>          If SIDE = 'L', IFAILR is not referenced.
00208 *> \endverbatim
00209 *>
00210 *> \param[out] INFO
00211 *> \verbatim
00212 *>          INFO is INTEGER
00213 *>          = 0:  successful exit
00214 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00215 *>          > 0:  if INFO = i, i is the number of eigenvectors which
00216 *>                failed to converge; see IFAILL and IFAILR for further
00217 *>                details.
00218 *> \endverbatim
00219 *
00220 *  Authors:
00221 *  ========
00222 *
00223 *> \author Univ. of Tennessee 
00224 *> \author Univ. of California Berkeley 
00225 *> \author Univ. of Colorado Denver 
00226 *> \author NAG Ltd. 
00227 *
00228 *> \date November 2011
00229 *
00230 *> \ingroup complex16OTHERcomputational
00231 *
00232 *> \par Further Details:
00233 *  =====================
00234 *>
00235 *> \verbatim
00236 *>
00237 *>  Each eigenvector is normalized so that the element of largest
00238 *>  magnitude has magnitude 1; here the magnitude of a complex number
00239 *>  (x,y) is taken to be |x|+|y|.
00240 *> \endverbatim
00241 *>
00242 *  =====================================================================
00243       SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
00244      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
00245      $                   IFAILR, INFO )
00246 *
00247 *  -- LAPACK computational routine (version 3.4.0) --
00248 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00249 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00250 *     November 2011
00251 *
00252 *     .. Scalar Arguments ..
00253       CHARACTER          EIGSRC, INITV, SIDE
00254       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
00255 *     ..
00256 *     .. Array Arguments ..
00257       LOGICAL            SELECT( * )
00258       INTEGER            IFAILL( * ), IFAILR( * )
00259       DOUBLE PRECISION   RWORK( * )
00260       COMPLEX*16         H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00261      $                   W( * ), WORK( * )
00262 *     ..
00263 *
00264 *  =====================================================================
00265 *
00266 *     .. Parameters ..
00267       COMPLEX*16         ZERO
00268       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
00269       DOUBLE PRECISION   RZERO
00270       PARAMETER          ( RZERO = 0.0D+0 )
00271 *     ..
00272 *     .. Local Scalars ..
00273       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
00274       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
00275       DOUBLE PRECISION   EPS3, HNORM, SMLNUM, ULP, UNFL
00276       COMPLEX*16         CDUM, WK
00277 *     ..
00278 *     .. External Functions ..
00279       LOGICAL            LSAME
00280       DOUBLE PRECISION   DLAMCH, ZLANHS
00281       EXTERNAL           LSAME, DLAMCH, ZLANHS
00282 *     ..
00283 *     .. External Subroutines ..
00284       EXTERNAL           XERBLA, ZLAEIN
00285 *     ..
00286 *     .. Intrinsic Functions ..
00287       INTRINSIC          ABS, DBLE, DIMAG, MAX
00288 *     ..
00289 *     .. Statement Functions ..
00290       DOUBLE PRECISION   CABS1
00291 *     ..
00292 *     .. Statement Function definitions ..
00293       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
00294 *     ..
00295 *     .. Executable Statements ..
00296 *
00297 *     Decode and test the input parameters.
00298 *
00299       BOTHV = LSAME( SIDE, 'B' )
00300       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
00301       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
00302 *
00303       FROMQR = LSAME( EIGSRC, 'Q' )
00304 *
00305       NOINIT = LSAME( INITV, 'N' )
00306 *
00307 *     Set M to the number of columns required to store the selected
00308 *     eigenvectors.
00309 *
00310       M = 0
00311       DO 10 K = 1, N
00312          IF( SELECT( K ) )
00313      $      M = M + 1
00314    10 CONTINUE
00315 *
00316       INFO = 0
00317       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
00318          INFO = -1
00319       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
00320          INFO = -2
00321       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
00322          INFO = -3
00323       ELSE IF( N.LT.0 ) THEN
00324          INFO = -5
00325       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
00326          INFO = -7
00327       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
00328          INFO = -10
00329       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
00330          INFO = -12
00331       ELSE IF( MM.LT.M ) THEN
00332          INFO = -13
00333       END IF
00334       IF( INFO.NE.0 ) THEN
00335          CALL XERBLA( 'ZHSEIN', -INFO )
00336          RETURN
00337       END IF
00338 *
00339 *     Quick return if possible.
00340 *
00341       IF( N.EQ.0 )
00342      $   RETURN
00343 *
00344 *     Set machine-dependent constants.
00345 *
00346       UNFL = DLAMCH( 'Safe minimum' )
00347       ULP = DLAMCH( 'Precision' )
00348       SMLNUM = UNFL*( N / ULP )
00349 *
00350       LDWORK = N
00351 *
00352       KL = 1
00353       KLN = 0
00354       IF( FROMQR ) THEN
00355          KR = 0
00356       ELSE
00357          KR = N
00358       END IF
00359       KS = 1
00360 *
00361       DO 100 K = 1, N
00362          IF( SELECT( K ) ) THEN
00363 *
00364 *           Compute eigenvector(s) corresponding to W(K).
00365 *
00366             IF( FROMQR ) THEN
00367 *
00368 *              If affiliation of eigenvalues is known, check whether
00369 *              the matrix splits.
00370 *
00371 *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
00372 *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
00373 *              KR = N).
00374 *
00375 *              Then inverse iteration can be performed with the
00376 *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
00377 *              the submatrix H(1:KR,1:KR) for a right eigenvector.
00378 *
00379                DO 20 I = K, KL + 1, -1
00380                   IF( H( I, I-1 ).EQ.ZERO )
00381      $               GO TO 30
00382    20          CONTINUE
00383    30          CONTINUE
00384                KL = I
00385                IF( K.GT.KR ) THEN
00386                   DO 40 I = K, N - 1
00387                      IF( H( I+1, I ).EQ.ZERO )
00388      $                  GO TO 50
00389    40             CONTINUE
00390    50             CONTINUE
00391                   KR = I
00392                END IF
00393             END IF
00394 *
00395             IF( KL.NE.KLN ) THEN
00396                KLN = KL
00397 *
00398 *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
00399 *              has not ben computed before.
00400 *
00401                HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
00402                IF( HNORM.GT.RZERO ) THEN
00403                   EPS3 = HNORM*ULP
00404                ELSE
00405                   EPS3 = SMLNUM
00406                END IF
00407             END IF
00408 *
00409 *           Perturb eigenvalue if it is close to any previous
00410 *           selected eigenvalues affiliated to the submatrix
00411 *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
00412 *
00413             WK = W( K )
00414    60       CONTINUE
00415             DO 70 I = K - 1, KL, -1
00416                IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
00417                   WK = WK + EPS3
00418                   GO TO 60
00419                END IF
00420    70       CONTINUE
00421             W( K ) = WK
00422 *
00423             IF( LEFTV ) THEN
00424 *
00425 *              Compute left eigenvector.
00426 *
00427                CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
00428      $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
00429      $                      SMLNUM, IINFO )
00430                IF( IINFO.GT.0 ) THEN
00431                   INFO = INFO + 1
00432                   IFAILL( KS ) = K
00433                ELSE
00434                   IFAILL( KS ) = 0
00435                END IF
00436                DO 80 I = 1, KL - 1
00437                   VL( I, KS ) = ZERO
00438    80          CONTINUE
00439             END IF
00440             IF( RIGHTV ) THEN
00441 *
00442 *              Compute right eigenvector.
00443 *
00444                CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
00445      $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
00446                IF( IINFO.GT.0 ) THEN
00447                   INFO = INFO + 1
00448                   IFAILR( KS ) = K
00449                ELSE
00450                   IFAILR( KS ) = 0
00451                END IF
00452                DO 90 I = KR + 1, N
00453                   VR( I, KS ) = ZERO
00454    90          CONTINUE
00455             END IF
00456             KS = KS + 1
00457          END IF
00458   100 CONTINUE
00459 *
00460       RETURN
00461 *
00462 *     End of ZHSEIN
00463 *
00464       END
 All Files Functions