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