1*> \brief \b CTREVC
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTREVC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrevc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrevc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrevc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22*                          LDVR, MM, M, WORK, RWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          HOWMNY, SIDE
26*       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
27*       ..
28*       .. Array Arguments ..
29*       LOGICAL            SELECT( * )
30*       REAL               RWORK( * )
31*       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
32*      $                   WORK( * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> CTREVC computes some or all of the right and/or left eigenvectors of
42*> a complex upper triangular matrix T.
43*> Matrices of this type are produced by the Schur factorization of
44*> a complex general matrix:  A = Q*T*Q**H, as computed by CHSEQR.
45*>
46*> The right eigenvector x and the left eigenvector y of T corresponding
47*> to an eigenvalue w are defined by:
48*>
49*>              T*x = w*x,     (y**H)*T = w*(y**H)
50*>
51*> where y**H denotes the conjugate transpose of the vector y.
52*> The eigenvalues are not input to this routine, but are read directly
53*> from the diagonal of T.
54*>
55*> This routine returns the matrices X and/or Y of right and left
56*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
57*> input matrix.  If Q is the unitary factor that reduces a matrix A to
58*> Schur form T, then Q*X and Q*Y are the matrices of right and left
59*> eigenvectors of A.
60*> \endverbatim
61*
62*  Arguments:
63*  ==========
64*
65*> \param[in] SIDE
66*> \verbatim
67*>          SIDE is CHARACTER*1
68*>          = 'R':  compute right eigenvectors only;
69*>          = 'L':  compute left eigenvectors only;
70*>          = 'B':  compute both right and left eigenvectors.
71*> \endverbatim
72*>
73*> \param[in] HOWMNY
74*> \verbatim
75*>          HOWMNY is CHARACTER*1
76*>          = 'A':  compute all right and/or left eigenvectors;
77*>          = 'B':  compute all right and/or left eigenvectors,
78*>                  backtransformed using the matrices supplied in
79*>                  VR and/or VL;
80*>          = 'S':  compute selected right and/or left eigenvectors,
81*>                  as indicated by the logical array SELECT.
82*> \endverbatim
83*>
84*> \param[in] SELECT
85*> \verbatim
86*>          SELECT is LOGICAL array, dimension (N)
87*>          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
88*>          computed.
89*>          The eigenvector corresponding to the j-th eigenvalue is
90*>          computed if SELECT(j) = .TRUE..
91*>          Not referenced if HOWMNY = 'A' or 'B'.
92*> \endverbatim
93*>
94*> \param[in] N
95*> \verbatim
96*>          N is INTEGER
97*>          The order of the matrix T. N >= 0.
98*> \endverbatim
99*>
100*> \param[in,out] T
101*> \verbatim
102*>          T is COMPLEX array, dimension (LDT,N)
103*>          The upper triangular matrix T.  T is modified, but restored
104*>          on exit.
105*> \endverbatim
106*>
107*> \param[in] LDT
108*> \verbatim
109*>          LDT is INTEGER
110*>          The leading dimension of the array T. LDT >= max(1,N).
111*> \endverbatim
112*>
113*> \param[in,out] VL
114*> \verbatim
115*>          VL is COMPLEX array, dimension (LDVL,MM)
116*>          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
117*>          contain an N-by-N matrix Q (usually the unitary matrix Q of
118*>          Schur vectors returned by CHSEQR).
119*>          On exit, if SIDE = 'L' or 'B', VL contains:
120*>          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
121*>          if HOWMNY = 'B', the matrix Q*Y;
122*>          if HOWMNY = 'S', the left eigenvectors of T specified by
123*>                           SELECT, stored consecutively in the columns
124*>                           of VL, in the same order as their
125*>                           eigenvalues.
126*>          Not referenced if SIDE = 'R'.
127*> \endverbatim
128*>
129*> \param[in] LDVL
130*> \verbatim
131*>          LDVL is INTEGER
132*>          The leading dimension of the array VL.  LDVL >= 1, and if
133*>          SIDE = 'L' or 'B', LDVL >= N.
134*> \endverbatim
135*>
136*> \param[in,out] VR
137*> \verbatim
138*>          VR is COMPLEX array, dimension (LDVR,MM)
139*>          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
140*>          contain an N-by-N matrix Q (usually the unitary matrix Q of
141*>          Schur vectors returned by CHSEQR).
142*>          On exit, if SIDE = 'R' or 'B', VR contains:
143*>          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
144*>          if HOWMNY = 'B', the matrix Q*X;
145*>          if HOWMNY = 'S', the right eigenvectors of T specified by
146*>                           SELECT, stored consecutively in the columns
147*>                           of VR, in the same order as their
148*>                           eigenvalues.
149*>          Not referenced if SIDE = 'L'.
150*> \endverbatim
151*>
152*> \param[in] LDVR
153*> \verbatim
154*>          LDVR is INTEGER
155*>          The leading dimension of the array VR.  LDVR >= 1, and if
156*>          SIDE = 'R' or 'B'; LDVR >= N.
157*> \endverbatim
158*>
159*> \param[in] MM
160*> \verbatim
161*>          MM is INTEGER
162*>          The number of columns in the arrays VL and/or VR. MM >= M.
163*> \endverbatim
164*>
165*> \param[out] M
166*> \verbatim
167*>          M is INTEGER
168*>          The number of columns in the arrays VL and/or VR actually
169*>          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
170*>          is set to N.  Each selected eigenvector occupies one
171*>          column.
172*> \endverbatim
173*>
174*> \param[out] WORK
175*> \verbatim
176*>          WORK is COMPLEX array, dimension (2*N)
177*> \endverbatim
178*>
179*> \param[out] RWORK
180*> \verbatim
181*>          RWORK is REAL array, dimension (N)
182*> \endverbatim
183*>
184*> \param[out] INFO
185*> \verbatim
186*>          INFO is INTEGER
187*>          = 0:  successful exit
188*>          < 0:  if INFO = -i, the i-th argument had an illegal value
189*> \endverbatim
190*
191*  Authors:
192*  ========
193*
194*> \author Univ. of Tennessee
195*> \author Univ. of California Berkeley
196*> \author Univ. of Colorado Denver
197*> \author NAG Ltd.
198*
199*> \date November 2011
200*
201*> \ingroup complexOTHERcomputational
202*
203*> \par Further Details:
204*  =====================
205*>
206*> \verbatim
207*>
208*>  The algorithm used in this program is basically backward (forward)
209*>  substitution, with scaling to make the the code robust against
210*>  possible overflow.
211*>
212*>  Each eigenvector is normalized so that the element of largest
213*>  magnitude has magnitude 1; here the magnitude of a complex number
214*>  (x,y) is taken to be |x| + |y|.
215*> \endverbatim
216*>
217*  =====================================================================
218      SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
219     $                   LDVR, MM, M, WORK, RWORK, INFO )
220*
221*  -- LAPACK computational routine (version 3.4.0) --
222*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
223*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
224*     November 2011
225*
226*     .. Scalar Arguments ..
227      CHARACTER          HOWMNY, SIDE
228      INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
229*     ..
230*     .. Array Arguments ..
231      LOGICAL            SELECT( * )
232      REAL               RWORK( * )
233      COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
234     $                   WORK( * )
235*     ..
236*
237*  =====================================================================
238*
239*     .. Parameters ..
240      REAL               ZERO, ONE
241      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
242      COMPLEX            CMZERO, CMONE
243      PARAMETER          ( CMZERO = ( 0.0E+0, 0.0E+0 ),
244     $                   CMONE = ( 1.0E+0, 0.0E+0 ) )
245*     ..
246*     .. Local Scalars ..
247      LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
248      INTEGER            I, II, IS, J, K, KI
249      REAL               OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
250      COMPLEX            CDUM
251*     ..
252*     .. External Functions ..
253      LOGICAL            LSAME
254      INTEGER            ICAMAX
255      REAL               SCASUM, SLAMCH
256      EXTERNAL           LSAME, ICAMAX, SCASUM, SLAMCH
257*     ..
258*     .. External Subroutines ..
259      EXTERNAL           CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA
260*     ..
261*     .. Intrinsic Functions ..
262      INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL
263*     ..
264*     .. Statement Functions ..
265      REAL               CABS1
266*     ..
267*     .. Statement Function definitions ..
268      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
269*     ..
270*     .. Executable Statements ..
271*
272*     Decode and test the input parameters
273*
274      BOTHV = LSAME( SIDE, 'B' )
275      RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
276      LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
277*
278      ALLV = LSAME( HOWMNY, 'A' )
279      OVER = LSAME( HOWMNY, 'B' )
280      SOMEV = LSAME( HOWMNY, 'S' )
281*
282*     Set M to the number of columns required to store the selected
283*     eigenvectors.
284*
285      IF( SOMEV ) THEN
286         M = 0
287         DO 10 J = 1, N
288            IF( SELECT( J ) )
289     $         M = M + 1
290   10    CONTINUE
291      ELSE
292         M = N
293      END IF
294*
295      INFO = 0
296      IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
297         INFO = -1
298      ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
299         INFO = -2
300      ELSE IF( N.LT.0 ) THEN
301         INFO = -4
302      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
303         INFO = -6
304      ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
305         INFO = -8
306      ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
307         INFO = -10
308      ELSE IF( MM.LT.M ) THEN
309         INFO = -11
310      END IF
311      IF( INFO.NE.0 ) THEN
312         CALL XERBLA( 'CTREVC', -INFO )
313         RETURN
314      END IF
315*
316*     Quick return if possible.
317*
318      IF( N.EQ.0 )
319     $   RETURN
320*
321*     Set the constants to control overflow.
322*
323      UNFL = SLAMCH( 'Safe minimum' )
324      OVFL = ONE / UNFL
325      CALL SLABAD( UNFL, OVFL )
326      ULP = SLAMCH( 'Precision' )
327      SMLNUM = UNFL*( N / ULP )
328*
329*     Store the diagonal elements of T in working array WORK.
330*
331      DO 20 I = 1, N
332         WORK( I+N ) = T( I, I )
333   20 CONTINUE
334*
335*     Compute 1-norm of each column of strictly upper triangular
336*     part of T to control overflow in triangular solver.
337*
338      RWORK( 1 ) = ZERO
339      DO 30 J = 2, N
340         RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 )
341   30 CONTINUE
342*
343      IF( RIGHTV ) THEN
344*
345*        Compute right eigenvectors.
346*
347         IS = M
348         DO 80 KI = N, 1, -1
349*
350            IF( SOMEV ) THEN
351               IF( .NOT.SELECT( KI ) )
352     $            GO TO 80
353            END IF
354            SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
355*
356            WORK( 1 ) = CMONE
357*
358*           Form right-hand side.
359*
360            DO 40 K = 1, KI - 1
361               WORK( K ) = -T( K, KI )
362   40       CONTINUE
363*
364*           Solve the triangular system:
365*              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
366*
367            DO 50 K = 1, KI - 1
368               T( K, K ) = T( K, K ) - T( KI, KI )
369               IF( CABS1( T( K, K ) ).LT.SMIN )
370     $            T( K, K ) = SMIN
371   50       CONTINUE
372*
373            IF( KI.GT.1 ) THEN
374               CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
375     $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
376     $                      INFO )
377               WORK( KI ) = SCALE
378            END IF
379*
380*           Copy the vector x or Q*x to VR and normalize.
381*
382            IF( .NOT.OVER ) THEN
383               CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
384*
385               II = ICAMAX( KI, VR( 1, IS ), 1 )
386               REMAX = ONE / CABS1( VR( II, IS ) )
387               CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 )
388*
389               DO 60 K = KI + 1, N
390                  VR( K, IS ) = CMZERO
391   60          CONTINUE
392            ELSE
393               IF( KI.GT.1 )
394     $            CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
395     $                        1, CMPLX( SCALE ), VR( 1, KI ), 1 )
396*
397               II = ICAMAX( N, VR( 1, KI ), 1 )
398               REMAX = ONE / CABS1( VR( II, KI ) )
399               CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 )
400            END IF
401*
402*           Set back the original diagonal elements of T.
403*
404            DO 70 K = 1, KI - 1
405               T( K, K ) = WORK( K+N )
406   70       CONTINUE
407*
408            IS = IS - 1
409   80    CONTINUE
410      END IF
411*
412      IF( LEFTV ) THEN
413*
414*        Compute left eigenvectors.
415*
416         IS = 1
417         DO 130 KI = 1, N
418*
419            IF( SOMEV ) THEN
420               IF( .NOT.SELECT( KI ) )
421     $            GO TO 130
422            END IF
423            SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
424*
425            WORK( N ) = CMONE
426*
427*           Form right-hand side.
428*
429            DO 90 K = KI + 1, N
430               WORK( K ) = -CONJG( T( KI, K ) )
431   90       CONTINUE
432*
433*           Solve the triangular system:
434*              (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK.
435*
436            DO 100 K = KI + 1, N
437               T( K, K ) = T( K, K ) - T( KI, KI )
438               IF( CABS1( T( K, K ) ).LT.SMIN )
439     $            T( K, K ) = SMIN
440  100       CONTINUE
441*
442            IF( KI.LT.N ) THEN
443               CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
444     $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
445     $                      WORK( KI+1 ), SCALE, RWORK, INFO )
446               WORK( KI ) = SCALE
447            END IF
448*
449*           Copy the vector x or Q*x to VL and normalize.
450*
451            IF( .NOT.OVER ) THEN
452               CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
453*
454               II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
455               REMAX = ONE / CABS1( VL( II, IS ) )
456               CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
457*
458               DO 110 K = 1, KI - 1
459                  VL( K, IS ) = CMZERO
460  110          CONTINUE
461            ELSE
462               IF( KI.LT.N )
463     $            CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
464     $                        WORK( KI+1 ), 1, CMPLX( SCALE ),
465     $                        VL( 1, KI ), 1 )
466*
467               II = ICAMAX( N, VL( 1, KI ), 1 )
468               REMAX = ONE / CABS1( VL( II, KI ) )
469               CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 )
470            END IF
471*
472*           Set back the original diagonal elements of T.
473*
474            DO 120 K = KI + 1, N
475               T( K, K ) = WORK( K+N )
476  120       CONTINUE
477*
478            IS = IS + 1
479  130    CONTINUE
480      END IF
481*
482      RETURN
483*
484*     End of CTREVC
485*
486      END
487c $Id$
488