LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
chsein.f
Go to the documentation of this file.
00001 *> \brief \b CHSEIN
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CHSEIN + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chsein.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chsein.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chsein.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CHSEIN( 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 *       REAL               RWORK( * )
00033 *       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00034 *      $                   W( * ), WORK( * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N)
00183 *> \endverbatim
00184 *>
00185 *> \param[out] RWORK
00186 *> \verbatim
00187 *>          RWORK is REAL 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 complexOTHERcomputational
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 CHSEIN( 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       REAL               RWORK( * )
00260       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00261      $                   W( * ), WORK( * )
00262 *     ..
00263 *
00264 *  =====================================================================
00265 *
00266 *     .. Parameters ..
00267       COMPLEX            ZERO
00268       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
00269       REAL               RZERO
00270       PARAMETER          ( RZERO = 0.0E+0 )
00271 *     ..
00272 *     .. Local Scalars ..
00273       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
00274       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
00275       REAL               EPS3, HNORM, SMLNUM, ULP, UNFL
00276       COMPLEX            CDUM, WK
00277 *     ..
00278 *     .. External Functions ..
00279       LOGICAL            LSAME
00280       REAL               CLANHS, SLAMCH
00281       EXTERNAL           LSAME, CLANHS, SLAMCH
00282 *     ..
00283 *     .. External Subroutines ..
00284       EXTERNAL           CLAEIN, XERBLA
00285 *     ..
00286 *     .. Intrinsic Functions ..
00287       INTRINSIC          ABS, AIMAG, MAX, REAL
00288 *     ..
00289 *     .. Statement Functions ..
00290       REAL               CABS1
00291 *     ..
00292 *     .. Statement Function definitions ..
00293       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( 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( 'CHSEIN', -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 = SLAMCH( 'Safe minimum' )
00347       ULP = SLAMCH( '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 = CLANHS( '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 CLAEIN( .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 CLAEIN( .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 CHSEIN
00463 *
00464       END
 All Files Functions