001:       SUBROUTINE CHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
002:      $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
003:      $                   IFAILR, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          EIGSRC, INITV, SIDE
012:       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
013: *     ..
014: *     .. Array Arguments ..
015:       LOGICAL            SELECT( * )
016:       INTEGER            IFAILL( * ), IFAILR( * )
017:       REAL               RWORK( * )
018:       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
019:      $                   W( * ), WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  CHSEIN uses inverse iteration to find specified right and/or left
026: *  eigenvectors of a complex upper Hessenberg matrix H.
027: *
028: *  The right eigenvector x and the left eigenvector y of the matrix H
029: *  corresponding to an eigenvalue w are defined by:
030: *
031: *               H * x = w * x,     y**h * H = w * y**h
032: *
033: *  where y**h denotes the conjugate transpose of the vector y.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  SIDE    (input) CHARACTER*1
039: *          = 'R': compute right eigenvectors only;
040: *          = 'L': compute left eigenvectors only;
041: *          = 'B': compute both right and left eigenvectors.
042: *
043: *  EIGSRC  (input) CHARACTER*1
044: *          Specifies the source of eigenvalues supplied in W:
045: *          = 'Q': the eigenvalues were found using CHSEQR; thus, if
046: *                 H has zero subdiagonal elements, and so is
047: *                 block-triangular, then the j-th eigenvalue can be
048: *                 assumed to be an eigenvalue of the block containing
049: *                 the j-th row/column.  This property allows CHSEIN to
050: *                 perform inverse iteration on just one diagonal block.
051: *          = 'N': no assumptions are made on the correspondence
052: *                 between eigenvalues and diagonal blocks.  In this
053: *                 case, CHSEIN must always perform inverse iteration
054: *                 using the whole matrix H.
055: *
056: *  INITV   (input) CHARACTER*1
057: *          = 'N': no initial vectors are supplied;
058: *          = 'U': user-supplied initial vectors are stored in the arrays
059: *                 VL and/or VR.
060: *
061: *  SELECT  (input) LOGICAL array, dimension (N)
062: *          Specifies the eigenvectors to be computed. To select the
063: *          eigenvector corresponding to the eigenvalue W(j),
064: *          SELECT(j) must be set to .TRUE..
065: *
066: *  N       (input) INTEGER
067: *          The order of the matrix H.  N >= 0.
068: *
069: *  H       (input) COMPLEX array, dimension (LDH,N)
070: *          The upper Hessenberg matrix H.
071: *
072: *  LDH     (input) INTEGER
073: *          The leading dimension of the array H.  LDH >= max(1,N).
074: *
075: *  W       (input/output) COMPLEX array, dimension (N)
076: *          On entry, the eigenvalues of H.
077: *          On exit, the real parts of W may have been altered since
078: *          close eigenvalues are perturbed slightly in searching for
079: *          independent eigenvectors.
080: *
081: *  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
082: *          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
083: *          contain starting vectors for the inverse iteration for the
084: *          left eigenvectors; the starting vector for each eigenvector
085: *          must be in the same column in which the eigenvector will be
086: *          stored.
087: *          On exit, if SIDE = 'L' or 'B', the left eigenvectors
088: *          specified by SELECT will be stored consecutively in the
089: *          columns of VL, in the same order as their eigenvalues.
090: *          If SIDE = 'R', VL is not referenced.
091: *
092: *  LDVL    (input) INTEGER
093: *          The leading dimension of the array VL.
094: *          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
095: *
096: *  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
097: *          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
098: *          contain starting vectors for the inverse iteration for the
099: *          right eigenvectors; the starting vector for each eigenvector
100: *          must be in the same column in which the eigenvector will be
101: *          stored.
102: *          On exit, if SIDE = 'R' or 'B', the right eigenvectors
103: *          specified by SELECT will be stored consecutively in the
104: *          columns of VR, in the same order as their eigenvalues.
105: *          If SIDE = 'L', VR is not referenced.
106: *
107: *  LDVR    (input) INTEGER
108: *          The leading dimension of the array VR.
109: *          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
110: *
111: *  MM      (input) INTEGER
112: *          The number of columns in the arrays VL and/or VR. MM >= M.
113: *
114: *  M       (output) INTEGER
115: *          The number of columns in the arrays VL and/or VR required to
116: *          store the eigenvectors (= the number of .TRUE. elements in
117: *          SELECT).
118: *
119: *  WORK    (workspace) COMPLEX array, dimension (N*N)
120: *
121: *  RWORK   (workspace) REAL array, dimension (N)
122: *
123: *  IFAILL  (output) INTEGER array, dimension (MM)
124: *          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
125: *          eigenvector in the i-th column of VL (corresponding to the
126: *          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
127: *          eigenvector converged satisfactorily.
128: *          If SIDE = 'R', IFAILL is not referenced.
129: *
130: *  IFAILR  (output) INTEGER array, dimension (MM)
131: *          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
132: *          eigenvector in the i-th column of VR (corresponding to the
133: *          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
134: *          eigenvector converged satisfactorily.
135: *          If SIDE = 'L', IFAILR is not referenced.
136: *
137: *  INFO    (output) INTEGER
138: *          = 0:  successful exit
139: *          < 0:  if INFO = -i, the i-th argument had an illegal value
140: *          > 0:  if INFO = i, i is the number of eigenvectors which
141: *                failed to converge; see IFAILL and IFAILR for further
142: *                details.
143: *
144: *  Further Details
145: *  ===============
146: *
147: *  Each eigenvector is normalized so that the element of largest
148: *  magnitude has magnitude 1; here the magnitude of a complex number
149: *  (x,y) is taken to be |x|+|y|.
150: *
151: *  =====================================================================
152: *
153: *     .. Parameters ..
154:       COMPLEX            ZERO
155:       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
156:       REAL               RZERO
157:       PARAMETER          ( RZERO = 0.0E+0 )
158: *     ..
159: *     .. Local Scalars ..
160:       LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
161:       INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
162:       REAL               EPS3, HNORM, SMLNUM, ULP, UNFL
163:       COMPLEX            CDUM, WK
164: *     ..
165: *     .. External Functions ..
166:       LOGICAL            LSAME
167:       REAL               CLANHS, SLAMCH
168:       EXTERNAL           LSAME, CLANHS, SLAMCH
169: *     ..
170: *     .. External Subroutines ..
171:       EXTERNAL           CLAEIN, XERBLA
172: *     ..
173: *     .. Intrinsic Functions ..
174:       INTRINSIC          ABS, AIMAG, MAX, REAL
175: *     ..
176: *     .. Statement Functions ..
177:       REAL               CABS1
178: *     ..
179: *     .. Statement Function definitions ..
180:       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
181: *     ..
182: *     .. Executable Statements ..
183: *
184: *     Decode and test the input parameters.
185: *
186:       BOTHV = LSAME( SIDE, 'B' )
187:       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
188:       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
189: *
190:       FROMQR = LSAME( EIGSRC, 'Q' )
191: *
192:       NOINIT = LSAME( INITV, 'N' )
193: *
194: *     Set M to the number of columns required to store the selected
195: *     eigenvectors.
196: *
197:       M = 0
198:       DO 10 K = 1, N
199:          IF( SELECT( K ) )
200:      $      M = M + 1
201:    10 CONTINUE
202: *
203:       INFO = 0
204:       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
205:          INFO = -1
206:       ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
207:          INFO = -2
208:       ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
209:          INFO = -3
210:       ELSE IF( N.LT.0 ) THEN
211:          INFO = -5
212:       ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
213:          INFO = -7
214:       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
215:          INFO = -10
216:       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
217:          INFO = -12
218:       ELSE IF( MM.LT.M ) THEN
219:          INFO = -13
220:       END IF
221:       IF( INFO.NE.0 ) THEN
222:          CALL XERBLA( 'CHSEIN', -INFO )
223:          RETURN
224:       END IF
225: *
226: *     Quick return if possible.
227: *
228:       IF( N.EQ.0 )
229:      $   RETURN
230: *
231: *     Set machine-dependent constants.
232: *
233:       UNFL = SLAMCH( 'Safe minimum' )
234:       ULP = SLAMCH( 'Precision' )
235:       SMLNUM = UNFL*( N / ULP )
236: *
237:       LDWORK = N
238: *
239:       KL = 1
240:       KLN = 0
241:       IF( FROMQR ) THEN
242:          KR = 0
243:       ELSE
244:          KR = N
245:       END IF
246:       KS = 1
247: *
248:       DO 100 K = 1, N
249:          IF( SELECT( K ) ) THEN
250: *
251: *           Compute eigenvector(s) corresponding to W(K).
252: *
253:             IF( FROMQR ) THEN
254: *
255: *              If affiliation of eigenvalues is known, check whether
256: *              the matrix splits.
257: *
258: *              Determine KL and KR such that 1 <= KL <= K <= KR <= N
259: *              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
260: *              KR = N).
261: *
262: *              Then inverse iteration can be performed with the
263: *              submatrix H(KL:N,KL:N) for a left eigenvector, and with
264: *              the submatrix H(1:KR,1:KR) for a right eigenvector.
265: *
266:                DO 20 I = K, KL + 1, -1
267:                   IF( H( I, I-1 ).EQ.ZERO )
268:      $               GO TO 30
269:    20          CONTINUE
270:    30          CONTINUE
271:                KL = I
272:                IF( K.GT.KR ) THEN
273:                   DO 40 I = K, N - 1
274:                      IF( H( I+1, I ).EQ.ZERO )
275:      $                  GO TO 50
276:    40             CONTINUE
277:    50             CONTINUE
278:                   KR = I
279:                END IF
280:             END IF
281: *
282:             IF( KL.NE.KLN ) THEN
283:                KLN = KL
284: *
285: *              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
286: *              has not ben computed before.
287: *
288:                HNORM = CLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
289:                IF( HNORM.GT.RZERO ) THEN
290:                   EPS3 = HNORM*ULP
291:                ELSE
292:                   EPS3 = SMLNUM
293:                END IF
294:             END IF
295: *
296: *           Perturb eigenvalue if it is close to any previous
297: *           selected eigenvalues affiliated to the submatrix
298: *           H(KL:KR,KL:KR). Close roots are modified by EPS3.
299: *
300:             WK = W( K )
301:    60       CONTINUE
302:             DO 70 I = K - 1, KL, -1
303:                IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
304:                   WK = WK + EPS3
305:                   GO TO 60
306:                END IF
307:    70       CONTINUE
308:             W( K ) = WK
309: *
310:             IF( LEFTV ) THEN
311: *
312: *              Compute left eigenvector.
313: *
314:                CALL CLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
315:      $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
316:      $                      SMLNUM, IINFO )
317:                IF( IINFO.GT.0 ) THEN
318:                   INFO = INFO + 1
319:                   IFAILL( KS ) = K
320:                ELSE
321:                   IFAILL( KS ) = 0
322:                END IF
323:                DO 80 I = 1, KL - 1
324:                   VL( I, KS ) = ZERO
325:    80          CONTINUE
326:             END IF
327:             IF( RIGHTV ) THEN
328: *
329: *              Compute right eigenvector.
330: *
331:                CALL CLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
332:      $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
333:                IF( IINFO.GT.0 ) THEN
334:                   INFO = INFO + 1
335:                   IFAILR( KS ) = K
336:                ELSE
337:                   IFAILR( KS ) = 0
338:                END IF
339:                DO 90 I = KR + 1, N
340:                   VR( I, KS ) = ZERO
341:    90          CONTINUE
342:             END IF
343:             KS = KS + 1
344:          END IF
345:   100 CONTINUE
346: *
347:       RETURN
348: *
349: *     End of CHSEIN
350: *
351:       END
352: