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