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