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