1*> \brief \b CHSEIN
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CHSEIN + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chsein.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chsein.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chsein.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
22*                          LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
23*                          IFAILR, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          EIGSRC, INITV, SIDE
27*       INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
28*       ..
29*       .. Array Arguments ..
30*       LOGICAL            SELECT( * )
31*       INTEGER            IFAILL( * ), IFAILR( * )
32*       REAL               RWORK( * )
33*       COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
34*      $                   W( * ), WORK( * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> CHSEIN uses inverse iteration to find specified right and/or left
44*> eigenvectors of a complex upper Hessenberg matrix H.
45*>
46*> The right eigenvector x and the left eigenvector y of the matrix H
47*> corresponding to an eigenvalue w are defined by:
48*>
49*>              H * x = w * x,     y**h * H = w * y**h
50*>
51*> where y**h denotes the conjugate transpose of the vector y.
52*> \endverbatim
53*
54*  Arguments:
55*  ==========
56*
57*> \param[in] SIDE
58*> \verbatim
59*>          SIDE is CHARACTER*1
60*>          = 'R': compute right eigenvectors only;
61*>          = 'L': compute left eigenvectors only;
62*>          = 'B': compute both right and left eigenvectors.
63*> \endverbatim
64*>
65*> \param[in] EIGSRC
66*> \verbatim
67*>          EIGSRC is CHARACTER*1
68*>          Specifies the source of eigenvalues supplied in W:
69*>          = 'Q': the eigenvalues were found using CHSEQR; thus, if
70*>                 H has zero subdiagonal elements, and so is
71*>                 block-triangular, then the j-th eigenvalue can be
72*>                 assumed to be an eigenvalue of the block containing
73*>                 the j-th row/column.  This property allows CHSEIN to
74*>                 perform inverse iteration on just one diagonal block.
75*>          = 'N': no assumptions are made on the correspondence
76*>                 between eigenvalues and diagonal blocks.  In this
77*>                 case, CHSEIN must always perform inverse iteration
78*>                 using the whole matrix H.
79*> \endverbatim
80*>
81*> \param[in] INITV
82*> \verbatim
83*>          INITV is CHARACTER*1
84*>          = 'N': no initial vectors are supplied;
85*>          = 'U': user-supplied initial vectors are stored in the arrays
86*>                 VL and/or VR.
87*> \endverbatim
88*>
89*> \param[in] SELECT
90*> \verbatim
91*>          SELECT is LOGICAL array, dimension (N)
92*>          Specifies the eigenvectors to be computed. To select the
93*>          eigenvector corresponding to the eigenvalue W(j),
94*>          SELECT(j) must be set to .TRUE..
95*> \endverbatim
96*>
97*> \param[in] N
98*> \verbatim
99*>          N is INTEGER
100*>          The order of the matrix H.  N >= 0.
101*> \endverbatim
102*>
103*> \param[in] H
104*> \verbatim
105*>          H is COMPLEX array, dimension (LDH,N)
106*>          The upper Hessenberg matrix H.
107*>          If a NaN is detected in H, the routine will return with INFO=-6.
108*> \endverbatim
109*>
110*> \param[in] LDH
111*> \verbatim
112*>          LDH is INTEGER
113*>          The leading dimension of the array H.  LDH >= max(1,N).
114*> \endverbatim
115*>
116*> \param[in,out] W
117*> \verbatim
118*>          W is COMPLEX array, dimension (N)
119*>          On entry, the eigenvalues of H.
120*>          On exit, the real parts of W may have been altered since
121*>          close eigenvalues are perturbed slightly in searching for
122*>          independent eigenvectors.
123*> \endverbatim
124*>
125*> \param[in,out] VL
126*> \verbatim
127*>          VL is COMPLEX array, dimension (LDVL,MM)
128*>          On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
129*>          contain starting vectors for the inverse iteration for the
130*>          left eigenvectors; the starting vector for each eigenvector
131*>          must be in the same column in which the eigenvector will be
132*>          stored.
133*>          On exit, if SIDE = 'L' or 'B', the left eigenvectors
134*>          specified by SELECT will be stored consecutively in the
135*>          columns of VL, in the same order as their eigenvalues.
136*>          If SIDE = 'R', VL is not referenced.
137*> \endverbatim
138*>
139*> \param[in] LDVL
140*> \verbatim
141*>          LDVL is INTEGER
142*>          The leading dimension of the array VL.
143*>          LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
144*> \endverbatim
145*>
146*> \param[in,out] VR
147*> \verbatim
148*>          VR is COMPLEX array, dimension (LDVR,MM)
149*>          On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
150*>          contain starting vectors for the inverse iteration for the
151*>          right eigenvectors; the starting vector for each eigenvector
152*>          must be in the same column in which the eigenvector will be
153*>          stored.
154*>          On exit, if SIDE = 'R' or 'B', the right eigenvectors
155*>          specified by SELECT will be stored consecutively in the
156*>          columns of VR, in the same order as their eigenvalues.
157*>          If SIDE = 'L', VR is not referenced.
158*> \endverbatim
159*>
160*> \param[in] LDVR
161*> \verbatim
162*>          LDVR is INTEGER
163*>          The leading dimension of the array VR.
164*>          LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
165*> \endverbatim
166*>
167*> \param[in] MM
168*> \verbatim
169*>          MM is INTEGER
170*>          The number of columns in the arrays VL and/or VR. MM >= M.
171*> \endverbatim
172*>
173*> \param[out] M
174*> \verbatim
175*>          M is INTEGER
176*>          The number of columns in the arrays VL and/or VR required to
177*>          store the eigenvectors (= the number of .TRUE. elements in
178*>          SELECT).
179*> \endverbatim
180*>
181*> \param[out] WORK
182*> \verbatim
183*>          WORK is COMPLEX array, dimension (N*N)
184*> \endverbatim
185*>
186*> \param[out] RWORK
187*> \verbatim
188*>          RWORK is REAL array, dimension (N)
189*> \endverbatim
190*>
191*> \param[out] IFAILL
192*> \verbatim
193*>          IFAILL is INTEGER array, dimension (MM)
194*>          If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
195*>          eigenvector in the i-th column of VL (corresponding to the
196*>          eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
197*>          eigenvector converged satisfactorily.
198*>          If SIDE = 'R', IFAILL is not referenced.
199*> \endverbatim
200*>
201*> \param[out] IFAILR
202*> \verbatim
203*>          IFAILR is INTEGER array, dimension (MM)
204*>          If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
205*>          eigenvector in the i-th column of VR (corresponding to the
206*>          eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
207*>          eigenvector converged satisfactorily.
208*>          If SIDE = 'L', IFAILR is not referenced.
209*> \endverbatim
210*>
211*> \param[out] INFO
212*> \verbatim
213*>          INFO is INTEGER
214*>          = 0:  successful exit
215*>          < 0:  if INFO = -i, the i-th argument had an illegal value
216*>          > 0:  if INFO = i, i is the number of eigenvectors which
217*>                failed to converge; see IFAILL and IFAILR for further
218*>                details.
219*> \endverbatim
220*
221*  Authors:
222*  ========
223*
224*> \author Univ. of Tennessee
225*> \author Univ. of California Berkeley
226*> \author Univ. of Colorado Denver
227*> \author NAG Ltd.
228*
229*> \ingroup complexOTHERcomputational
230*
231*> \par Further Details:
232*  =====================
233*>
234*> \verbatim
235*>
236*>  Each eigenvector is normalized so that the element of largest
237*>  magnitude has magnitude 1; here the magnitude of a complex number
238*>  (x,y) is taken to be |x|+|y|.
239*> \endverbatim
240*>
241*  =====================================================================
242      SUBROUTINE CHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
243     $                   LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
244     $                   IFAILR, INFO )
245*
246*  -- LAPACK computational routine --
247*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
248*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
249*
250*     .. Scalar Arguments ..
251      CHARACTER          EIGSRC, INITV, SIDE
252      INTEGER            INFO, LDH, LDVL, LDVR, M, MM, N
253*     ..
254*     .. Array Arguments ..
255      LOGICAL            SELECT( * )
256      INTEGER            IFAILL( * ), IFAILR( * )
257      REAL               RWORK( * )
258      COMPLEX            H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
259     $                   W( * ), WORK( * )
260*     ..
261*
262*  =====================================================================
263*
264*     .. Parameters ..
265      COMPLEX            ZERO
266      PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
267      REAL               RZERO
268      PARAMETER          ( RZERO = 0.0E+0 )
269*     ..
270*     .. Local Scalars ..
271      LOGICAL            BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
272      INTEGER            I, IINFO, K, KL, KLN, KR, KS, LDWORK
273      REAL               EPS3, HNORM, SMLNUM, ULP, UNFL
274      COMPLEX            CDUM, WK
275*     ..
276*     .. External Functions ..
277      LOGICAL            LSAME, SISNAN
278      REAL               CLANHS, SLAMCH
279      EXTERNAL           LSAME, CLANHS, SLAMCH, SISNAN
280*     ..
281*     .. External Subroutines ..
282      EXTERNAL           CLAEIN, XERBLA
283*     ..
284*     .. Intrinsic Functions ..
285      INTRINSIC          ABS, AIMAG, MAX, REAL
286*     ..
287*     .. Statement Functions ..
288      REAL               CABS1
289*     ..
290*     .. Statement Function definitions ..
291      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
292*     ..
293*     .. Executable Statements ..
294*
295*     Decode and test the input parameters.
296*
297      BOTHV = LSAME( SIDE, 'B' )
298      RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
299      LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
300*
301      FROMQR = LSAME( EIGSRC, 'Q' )
302*
303      NOINIT = LSAME( INITV, 'N' )
304*
305*     Set M to the number of columns required to store the selected
306*     eigenvectors.
307*
308      M = 0
309      DO 10 K = 1, N
310         IF( SELECT( K ) )
311     $      M = M + 1
312   10 CONTINUE
313*
314      INFO = 0
315      IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
316         INFO = -1
317      ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
318         INFO = -2
319      ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
320         INFO = -3
321      ELSE IF( N.LT.0 ) THEN
322         INFO = -5
323      ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
324         INFO = -7
325      ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
326         INFO = -10
327      ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
328         INFO = -12
329      ELSE IF( MM.LT.M ) THEN
330         INFO = -13
331      END IF
332      IF( INFO.NE.0 ) THEN
333         CALL XERBLA( 'CHSEIN', -INFO )
334         RETURN
335      END IF
336*
337*     Quick return if possible.
338*
339      IF( N.EQ.0 )
340     $   RETURN
341*
342*     Set machine-dependent constants.
343*
344      UNFL = SLAMCH( 'Safe minimum' )
345      ULP = SLAMCH( 'Precision' )
346      SMLNUM = UNFL*( N / ULP )
347*
348      LDWORK = N
349*
350      KL = 1
351      KLN = 0
352      IF( FROMQR ) THEN
353         KR = 0
354      ELSE
355         KR = N
356      END IF
357      KS = 1
358*
359      DO 100 K = 1, N
360         IF( SELECT( K ) ) THEN
361*
362*           Compute eigenvector(s) corresponding to W(K).
363*
364            IF( FROMQR ) THEN
365*
366*              If affiliation of eigenvalues is known, check whether
367*              the matrix splits.
368*
369*              Determine KL and KR such that 1 <= KL <= K <= KR <= N
370*              and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
371*              KR = N).
372*
373*              Then inverse iteration can be performed with the
374*              submatrix H(KL:N,KL:N) for a left eigenvector, and with
375*              the submatrix H(1:KR,1:KR) for a right eigenvector.
376*
377               DO 20 I = K, KL + 1, -1
378                  IF( H( I, I-1 ).EQ.ZERO )
379     $               GO TO 30
380   20          CONTINUE
381   30          CONTINUE
382               KL = I
383               IF( K.GT.KR ) THEN
384                  DO 40 I = K, N - 1
385                     IF( H( I+1, I ).EQ.ZERO )
386     $                  GO TO 50
387   40             CONTINUE
388   50             CONTINUE
389                  KR = I
390               END IF
391            END IF
392*
393            IF( KL.NE.KLN ) THEN
394               KLN = KL
395*
396*              Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
397*              has not ben computed before.
398*
399               HNORM = CLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
400               IF( SISNAN( HNORM ) ) THEN
401                  INFO = -6
402                  RETURN
403               ELSE IF( (HNORM.GT.RZERO) ) THEN
404                  EPS3 = HNORM*ULP
405               ELSE
406                  EPS3 = SMLNUM
407               END IF
408            END IF
409*
410*           Perturb eigenvalue if it is close to any previous
411*           selected eigenvalues affiliated to the submatrix
412*           H(KL:KR,KL:KR). Close roots are modified by EPS3.
413*
414            WK = W( K )
415   60       CONTINUE
416            DO 70 I = K - 1, KL, -1
417               IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
418                  WK = WK + EPS3
419                  GO TO 60
420               END IF
421   70       CONTINUE
422            W( K ) = WK
423*
424            IF( LEFTV ) THEN
425*
426*              Compute left eigenvector.
427*
428               CALL CLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
429     $                      WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
430     $                      SMLNUM, IINFO )
431               IF( IINFO.GT.0 ) THEN
432                  INFO = INFO + 1
433                  IFAILL( KS ) = K
434               ELSE
435                  IFAILL( KS ) = 0
436               END IF
437               DO 80 I = 1, KL - 1
438                  VL( I, KS ) = ZERO
439   80          CONTINUE
440            END IF
441            IF( RIGHTV ) THEN
442*
443*              Compute right eigenvector.
444*
445               CALL CLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
446     $                      WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
447               IF( IINFO.GT.0 ) THEN
448                  INFO = INFO + 1
449                  IFAILR( KS ) = K
450               ELSE
451                  IFAILR( KS ) = 0
452               END IF
453               DO 90 I = KR + 1, N
454                  VR( I, KS ) = ZERO
455   90          CONTINUE
456            END IF
457            KS = KS + 1
458         END IF
459  100 CONTINUE
460*
461      RETURN
462*
463*     End of CHSEIN
464*
465      END
466