001:       SUBROUTINE CTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
002:      $                   LDVR, MM, M, WORK, RWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          HOWMNY, SIDE
010:       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
011: *     ..
012: *     .. Array Arguments ..
013:       LOGICAL            SELECT( * )
014:       REAL               RWORK( * )
015:       COMPLEX            T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
016:      $                   WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  CTREVC computes some or all of the right and/or left eigenvectors of
023: *  a complex upper triangular matrix T.
024: *  Matrices of this type are produced by the Schur factorization of
025: *  a complex general matrix:  A = Q*T*Q**H, as computed by CHSEQR.
026: *  
027: *  The right eigenvector x and the left eigenvector y of T corresponding
028: *  to an eigenvalue w are defined by:
029: *  
030: *               T*x = w*x,     (y**H)*T = w*(y**H)
031: *  
032: *  where y**H denotes the conjugate transpose of the vector y.
033: *  The eigenvalues are not input to this routine, but are read directly
034: *  from the diagonal of T.
035: *  
036: *  This routine returns the matrices X and/or Y of right and left
037: *  eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
038: *  input matrix.  If Q is the unitary factor that reduces a matrix A to
039: *  Schur form T, then Q*X and Q*Y are the matrices of right and left
040: *  eigenvectors of A.
041: *
042: *  Arguments
043: *  =========
044: *
045: *  SIDE    (input) CHARACTER*1
046: *          = 'R':  compute right eigenvectors only;
047: *          = 'L':  compute left eigenvectors only;
048: *          = 'B':  compute both right and left eigenvectors.
049: *
050: *  HOWMNY  (input) CHARACTER*1
051: *          = 'A':  compute all right and/or left eigenvectors;
052: *          = 'B':  compute all right and/or left eigenvectors,
053: *                  backtransformed using the matrices supplied in
054: *                  VR and/or VL;
055: *          = 'S':  compute selected right and/or left eigenvectors,
056: *                  as indicated by the logical array SELECT.
057: *
058: *  SELECT  (input) LOGICAL array, dimension (N)
059: *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
060: *          computed.
061: *          The eigenvector corresponding to the j-th eigenvalue is
062: *          computed if SELECT(j) = .TRUE..
063: *          Not referenced if HOWMNY = 'A' or 'B'.
064: *
065: *  N       (input) INTEGER
066: *          The order of the matrix T. N >= 0.
067: *
068: *  T       (input/output) COMPLEX array, dimension (LDT,N)
069: *          The upper triangular matrix T.  T is modified, but restored
070: *          on exit.
071: *
072: *  LDT     (input) INTEGER
073: *          The leading dimension of the array T. LDT >= max(1,N).
074: *
075: *  VL      (input/output) COMPLEX array, dimension (LDVL,MM)
076: *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
077: *          contain an N-by-N matrix Q (usually the unitary matrix Q of
078: *          Schur vectors returned by CHSEQR).
079: *          On exit, if SIDE = 'L' or 'B', VL contains:
080: *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
081: *          if HOWMNY = 'B', the matrix Q*Y;
082: *          if HOWMNY = 'S', the left eigenvectors of T specified by
083: *                           SELECT, stored consecutively in the columns
084: *                           of VL, in the same order as their
085: *                           eigenvalues.
086: *          Not referenced if SIDE = 'R'.
087: *
088: *  LDVL    (input) INTEGER
089: *          The leading dimension of the array VL.  LDVL >= 1, and if
090: *          SIDE = 'L' or 'B', LDVL >= N.
091: *
092: *  VR      (input/output) COMPLEX array, dimension (LDVR,MM)
093: *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
094: *          contain an N-by-N matrix Q (usually the unitary matrix Q of
095: *          Schur vectors returned by CHSEQR).
096: *          On exit, if SIDE = 'R' or 'B', VR contains:
097: *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
098: *          if HOWMNY = 'B', the matrix Q*X;
099: *          if HOWMNY = 'S', the right eigenvectors of T specified by
100: *                           SELECT, stored consecutively in the columns
101: *                           of VR, in the same order as their
102: *                           eigenvalues.
103: *          Not referenced if SIDE = 'L'.
104: *
105: *  LDVR    (input) INTEGER
106: *          The leading dimension of the array VR.  LDVR >= 1, and if
107: *          SIDE = 'R' or 'B'; LDVR >= N.
108: *
109: *  MM      (input) INTEGER
110: *          The number of columns in the arrays VL and/or VR. MM >= M.
111: *
112: *  M       (output) INTEGER
113: *          The number of columns in the arrays VL and/or VR actually
114: *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
115: *          is set to N.  Each selected eigenvector occupies one
116: *          column.
117: *
118: *  WORK    (workspace) COMPLEX array, dimension (2*N)
119: *
120: *  RWORK   (workspace) REAL array, dimension (N)
121: *
122: *  INFO    (output) INTEGER
123: *          = 0:  successful exit
124: *          < 0:  if INFO = -i, the i-th argument had an illegal value
125: *
126: *  Further Details
127: *  ===============
128: *
129: *  The algorithm used in this program is basically backward (forward)
130: *  substitution, with scaling to make the the code robust against
131: *  possible overflow.
132: *
133: *  Each eigenvector is normalized so that the element of largest
134: *  magnitude has magnitude 1; here the magnitude of a complex number
135: *  (x,y) is taken to be |x| + |y|.
136: *
137: *  =====================================================================
138: *
139: *     .. Parameters ..
140:       REAL               ZERO, ONE
141:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
142:       COMPLEX            CMZERO, CMONE
143:       PARAMETER          ( CMZERO = ( 0.0E+0, 0.0E+0 ),
144:      $                   CMONE = ( 1.0E+0, 0.0E+0 ) )
145: *     ..
146: *     .. Local Scalars ..
147:       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
148:       INTEGER            I, II, IS, J, K, KI
149:       REAL               OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
150:       COMPLEX            CDUM
151: *     ..
152: *     .. External Functions ..
153:       LOGICAL            LSAME
154:       INTEGER            ICAMAX
155:       REAL               SCASUM, SLAMCH
156:       EXTERNAL           LSAME, ICAMAX, SCASUM, SLAMCH
157: *     ..
158: *     .. External Subroutines ..
159:       EXTERNAL           CCOPY, CGEMV, CLATRS, CSSCAL, SLABAD, XERBLA
160: *     ..
161: *     .. Intrinsic Functions ..
162:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL
163: *     ..
164: *     .. Statement Functions ..
165:       REAL               CABS1
166: *     ..
167: *     .. Statement Function definitions ..
168:       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
169: *     ..
170: *     .. Executable Statements ..
171: *
172: *     Decode and test the input parameters
173: *
174:       BOTHV = LSAME( SIDE, 'B' )
175:       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
176:       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
177: *
178:       ALLV = LSAME( HOWMNY, 'A' )
179:       OVER = LSAME( HOWMNY, 'B' )
180:       SOMEV = LSAME( HOWMNY, 'S' )
181: *
182: *     Set M to the number of columns required to store the selected
183: *     eigenvectors.
184: *
185:       IF( SOMEV ) THEN
186:          M = 0
187:          DO 10 J = 1, N
188:             IF( SELECT( J ) )
189:      $         M = M + 1
190:    10    CONTINUE
191:       ELSE
192:          M = N
193:       END IF
194: *
195:       INFO = 0
196:       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
197:          INFO = -1
198:       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
199:          INFO = -2
200:       ELSE IF( N.LT.0 ) THEN
201:          INFO = -4
202:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
203:          INFO = -6
204:       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
205:          INFO = -8
206:       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
207:          INFO = -10
208:       ELSE IF( MM.LT.M ) THEN
209:          INFO = -11
210:       END IF
211:       IF( INFO.NE.0 ) THEN
212:          CALL XERBLA( 'CTREVC', -INFO )
213:          RETURN
214:       END IF
215: *
216: *     Quick return if possible.
217: *
218:       IF( N.EQ.0 )
219:      $   RETURN
220: *
221: *     Set the constants to control overflow.
222: *
223:       UNFL = SLAMCH( 'Safe minimum' )
224:       OVFL = ONE / UNFL
225:       CALL SLABAD( UNFL, OVFL )
226:       ULP = SLAMCH( 'Precision' )
227:       SMLNUM = UNFL*( N / ULP )
228: *
229: *     Store the diagonal elements of T in working array WORK.
230: *
231:       DO 20 I = 1, N
232:          WORK( I+N ) = T( I, I )
233:    20 CONTINUE
234: *
235: *     Compute 1-norm of each column of strictly upper triangular
236: *     part of T to control overflow in triangular solver.
237: *
238:       RWORK( 1 ) = ZERO
239:       DO 30 J = 2, N
240:          RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 )
241:    30 CONTINUE
242: *
243:       IF( RIGHTV ) THEN
244: *
245: *        Compute right eigenvectors.
246: *
247:          IS = M
248:          DO 80 KI = N, 1, -1
249: *
250:             IF( SOMEV ) THEN
251:                IF( .NOT.SELECT( KI ) )
252:      $            GO TO 80
253:             END IF
254:             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
255: *
256:             WORK( 1 ) = CMONE
257: *
258: *           Form right-hand side.
259: *
260:             DO 40 K = 1, KI - 1
261:                WORK( K ) = -T( K, KI )
262:    40       CONTINUE
263: *
264: *           Solve the triangular system:
265: *              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
266: *
267:             DO 50 K = 1, KI - 1
268:                T( K, K ) = T( K, K ) - T( KI, KI )
269:                IF( CABS1( T( K, K ) ).LT.SMIN )
270:      $            T( K, K ) = SMIN
271:    50       CONTINUE
272: *
273:             IF( KI.GT.1 ) THEN
274:                CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
275:      $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
276:      $                      INFO )
277:                WORK( KI ) = SCALE
278:             END IF
279: *
280: *           Copy the vector x or Q*x to VR and normalize.
281: *
282:             IF( .NOT.OVER ) THEN
283:                CALL CCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
284: *
285:                II = ICAMAX( KI, VR( 1, IS ), 1 )
286:                REMAX = ONE / CABS1( VR( II, IS ) )
287:                CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 )
288: *
289:                DO 60 K = KI + 1, N
290:                   VR( K, IS ) = CMZERO
291:    60          CONTINUE
292:             ELSE
293:                IF( KI.GT.1 )
294:      $            CALL CGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
295:      $                        1, CMPLX( SCALE ), VR( 1, KI ), 1 )
296: *
297:                II = ICAMAX( N, VR( 1, KI ), 1 )
298:                REMAX = ONE / CABS1( VR( II, KI ) )
299:                CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 )
300:             END IF
301: *
302: *           Set back the original diagonal elements of T.
303: *
304:             DO 70 K = 1, KI - 1
305:                T( K, K ) = WORK( K+N )
306:    70       CONTINUE
307: *
308:             IS = IS - 1
309:    80    CONTINUE
310:       END IF
311: *
312:       IF( LEFTV ) THEN
313: *
314: *        Compute left eigenvectors.
315: *
316:          IS = 1
317:          DO 130 KI = 1, N
318: *
319:             IF( SOMEV ) THEN
320:                IF( .NOT.SELECT( KI ) )
321:      $            GO TO 130
322:             END IF
323:             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
324: *
325:             WORK( N ) = CMONE
326: *
327: *           Form right-hand side.
328: *
329:             DO 90 K = KI + 1, N
330:                WORK( K ) = -CONJG( T( KI, K ) )
331:    90       CONTINUE
332: *
333: *           Solve the triangular system:
334: *              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK.
335: *
336:             DO 100 K = KI + 1, N
337:                T( K, K ) = T( K, K ) - T( KI, KI )
338:                IF( CABS1( T( K, K ) ).LT.SMIN )
339:      $            T( K, K ) = SMIN
340:   100       CONTINUE
341: *
342:             IF( KI.LT.N ) THEN
343:                CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
344:      $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
345:      $                      WORK( KI+1 ), SCALE, RWORK, INFO )
346:                WORK( KI ) = SCALE
347:             END IF
348: *
349: *           Copy the vector x or Q*x to VL and normalize.
350: *
351:             IF( .NOT.OVER ) THEN
352:                CALL CCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
353: *
354:                II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
355:                REMAX = ONE / CABS1( VL( II, IS ) )
356:                CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
357: *
358:                DO 110 K = 1, KI - 1
359:                   VL( K, IS ) = CMZERO
360:   110          CONTINUE
361:             ELSE
362:                IF( KI.LT.N )
363:      $            CALL CGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
364:      $                        WORK( KI+1 ), 1, CMPLX( SCALE ),
365:      $                        VL( 1, KI ), 1 )
366: *
367:                II = ICAMAX( N, VL( 1, KI ), 1 )
368:                REMAX = ONE / CABS1( VL( II, KI ) )
369:                CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 )
370:             END IF
371: *
372: *           Set back the original diagonal elements of T.
373: *
374:             DO 120 K = KI + 1, N
375:                T( K, K ) = WORK( K+N )
376:   120       CONTINUE
377: *
378:             IS = IS + 1
379:   130    CONTINUE
380:       END IF
381: *
382:       RETURN
383: *
384: *     End of CTREVC
385: *
386:       END
387: