LAPACK 3.3.0

chsein.f

Go to the documentation of this file.
00001       SUBROUTINE CHSEIN( 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       REAL               RWORK( * )
00018       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
00019      $                   W( * ), WORK( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N)
00120 *
00121 *  RWORK   (workspace) REAL 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            ZERO
00155       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
00156       REAL               RZERO
00157       PARAMETER          ( RZERO = 0.0E+0 )
00158 *     ..
00159 *     .. Local Scalars ..
00160       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
00161       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
00162       REAL               EPS3, HNORM, SMLNUM, ULP, UNFL
00163       COMPLEX            CDUM, WK
00164 *     ..
00165 *     .. External Functions ..
00166       LOGICAL            LSAME
00167       REAL               CLANHS, SLAMCH
00168       EXTERNAL           LSAME, CLANHS, SLAMCH
00169 *     ..
00170 *     .. External Subroutines ..
00171       EXTERNAL           CLAEIN, XERBLA
00172 *     ..
00173 *     .. Intrinsic Functions ..
00174       INTRINSIC          ABS, AIMAG, MAX, REAL
00175 *     ..
00176 *     .. Statement Functions ..
00177       REAL               CABS1
00178 *     ..
00179 *     .. Statement Function definitions ..
00180       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( 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( 'CHSEIN', -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 = SLAMCH( 'Safe minimum' )
00234       ULP = SLAMCH( '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 = CLANHS( '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 CLAEIN( .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 CLAEIN( .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 CHSEIN
00350 *
00351       END
 All Files Functions