 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ ctrevc()

 subroutine ctrevc ( character SIDE, character HOWMNY, logical, dimension( * ) SELECT, integer N, complex, dimension( ldt, * ) T, integer LDT, complex, dimension( ldvl, * ) VL, integer LDVL, complex, dimension( ldvr, * ) VR, integer LDVR, integer MM, integer M, complex, dimension( * ) WORK, real, dimension( * ) RWORK, integer INFO )

CTREVC

Purpose:
``` CTREVC computes some or all of the right and/or left eigenvectors of
a complex upper triangular matrix T.
Matrices of this type are produced by the Schur factorization of
a complex general matrix:  A = Q*T*Q**H, as computed by CHSEQR.

The right eigenvector x and the left eigenvector y of T corresponding
to an eigenvalue w are defined by:

T*x = w*x,     (y**H)*T = w*(y**H)

where y**H denotes the conjugate transpose of the vector y.
The eigenvalues are not input to this routine, but are read directly
from the diagonal of T.

This routine returns the matrices X and/or Y of right and left
eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
input matrix.  If Q is the unitary factor that reduces a matrix A to
Schur form T, then Q*X and Q*Y are the matrices of right and left
eigenvectors of A.```
Parameters
 [in] SIDE ``` SIDE is CHARACTER*1 = 'R': compute right eigenvectors only; = 'L': compute left eigenvectors only; = 'B': compute both right and left eigenvectors.``` [in] HOWMNY ``` HOWMNY is CHARACTER*1 = 'A': compute all right and/or left eigenvectors; = 'B': compute all right and/or left eigenvectors, backtransformed using the matrices supplied in VR and/or VL; = 'S': compute selected right and/or left eigenvectors, as indicated by the logical array SELECT.``` [in] SELECT ``` SELECT is LOGICAL array, dimension (N) If HOWMNY = 'S', SELECT specifies the eigenvectors to be computed. The eigenvector corresponding to the j-th eigenvalue is computed if SELECT(j) = .TRUE.. Not referenced if HOWMNY = 'A' or 'B'.``` [in] N ``` N is INTEGER The order of the matrix T. N >= 0.``` [in,out] T ``` T is COMPLEX array, dimension (LDT,N) The upper triangular matrix T. T is modified, but restored on exit.``` [in] LDT ``` LDT is INTEGER The leading dimension of the array T. LDT >= max(1,N).``` [in,out] VL ``` VL is COMPLEX array, dimension (LDVL,MM) On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must contain an N-by-N matrix Q (usually the unitary matrix Q of Schur vectors returned by CHSEQR). On exit, if SIDE = 'L' or 'B', VL contains: if HOWMNY = 'A', the matrix Y of left eigenvectors of T; if HOWMNY = 'B', the matrix Q*Y; if HOWMNY = 'S', the left eigenvectors of T specified by SELECT, stored consecutively in the columns of VL, in the same order as their eigenvalues. Not referenced if SIDE = 'R'.``` [in] LDVL ``` LDVL is INTEGER The leading dimension of the array VL. LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.``` [in,out] VR ``` VR is COMPLEX array, dimension (LDVR,MM) On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must contain an N-by-N matrix Q (usually the unitary matrix Q of Schur vectors returned by CHSEQR). On exit, if SIDE = 'R' or 'B', VR contains: if HOWMNY = 'A', the matrix X of right eigenvectors of T; if HOWMNY = 'B', the matrix Q*X; if HOWMNY = 'S', the right eigenvectors of T specified by SELECT, stored consecutively in the columns of VR, in the same order as their eigenvalues. Not referenced if SIDE = 'L'.``` [in] LDVR ``` LDVR is INTEGER The leading dimension of the array VR. LDVR >= 1, and if SIDE = 'R' or 'B'; LDVR >= N.``` [in] MM ``` MM is INTEGER The number of columns in the arrays VL and/or VR. MM >= M.``` [out] M ``` M is INTEGER The number of columns in the arrays VL and/or VR actually used to store the eigenvectors. If HOWMNY = 'A' or 'B', M is set to N. Each selected eigenvector occupies one column.``` [out] WORK ` WORK is COMPLEX array, dimension (2*N)` [out] RWORK ` RWORK is REAL array, dimension (N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Further Details:
```  The algorithm used in this program is basically backward (forward)
substitution, with scaling to make the the code robust against
possible overflow.

Each eigenvector is normalized so that the element of largest
magnitude has magnitude 1; here the magnitude of a complex number
(x,y) is taken to be |x| + |y|.```

Definition at line 216 of file ctrevc.f.

218*
219* -- LAPACK computational routine --
220* -- LAPACK is a software package provided by Univ. of Tennessee, --
221* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222*
223* .. Scalar Arguments ..
224 CHARACTER HOWMNY, SIDE
225 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
226* ..
227* .. Array Arguments ..
228 LOGICAL SELECT( * )
229 REAL RWORK( * )
230 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
231 \$ WORK( * )
232* ..
233*
234* =====================================================================
235*
236* .. Parameters ..
237 REAL ZERO, ONE
238 parameter( zero = 0.0e+0, one = 1.0e+0 )
239 COMPLEX CMZERO, CMONE
240 parameter( cmzero = ( 0.0e+0, 0.0e+0 ),
241 \$ cmone = ( 1.0e+0, 0.0e+0 ) )
242* ..
243* .. Local Scalars ..
244 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
245 INTEGER I, II, IS, J, K, KI
246 REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
247 COMPLEX CDUM
248* ..
249* .. External Functions ..
250 LOGICAL LSAME
251 INTEGER ICAMAX
252 REAL SCASUM, SLAMCH
253 EXTERNAL lsame, icamax, scasum, slamch
254* ..
255* .. External Subroutines ..
256 EXTERNAL ccopy, cgemv, clatrs, csscal, slabad, xerbla
257* ..
258* .. Intrinsic Functions ..
259 INTRINSIC abs, aimag, cmplx, conjg, max, real
260* ..
261* .. Statement Functions ..
262 REAL CABS1
263* ..
264* .. Statement Function definitions ..
265 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
266* ..
267* .. Executable Statements ..
268*
269* Decode and test the input parameters
270*
271 bothv = lsame( side, 'B' )
272 rightv = lsame( side, 'R' ) .OR. bothv
273 leftv = lsame( side, 'L' ) .OR. bothv
274*
275 allv = lsame( howmny, 'A' )
276 over = lsame( howmny, 'B' )
277 somev = lsame( howmny, 'S' )
278*
279* Set M to the number of columns required to store the selected
280* eigenvectors.
281*
282 IF( somev ) THEN
283 m = 0
284 DO 10 j = 1, n
285 IF( SELECT( j ) )
286 \$ m = m + 1
287 10 CONTINUE
288 ELSE
289 m = n
290 END IF
291*
292 info = 0
293 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
294 info = -1
295 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
296 info = -2
297 ELSE IF( n.LT.0 ) THEN
298 info = -4
299 ELSE IF( ldt.LT.max( 1, n ) ) THEN
300 info = -6
301 ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
302 info = -8
303 ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
304 info = -10
305 ELSE IF( mm.LT.m ) THEN
306 info = -11
307 END IF
308 IF( info.NE.0 ) THEN
309 CALL xerbla( 'CTREVC', -info )
310 RETURN
311 END IF
312*
313* Quick return if possible.
314*
315 IF( n.EQ.0 )
316 \$ RETURN
317*
318* Set the constants to control overflow.
319*
320 unfl = slamch( 'Safe minimum' )
321 ovfl = one / unfl
322 CALL slabad( unfl, ovfl )
323 ulp = slamch( 'Precision' )
324 smlnum = unfl*( n / ulp )
325*
326* Store the diagonal elements of T in working array WORK.
327*
328 DO 20 i = 1, n
329 work( i+n ) = t( i, i )
330 20 CONTINUE
331*
332* Compute 1-norm of each column of strictly upper triangular
333* part of T to control overflow in triangular solver.
334*
335 rwork( 1 ) = zero
336 DO 30 j = 2, n
337 rwork( j ) = scasum( j-1, t( 1, j ), 1 )
338 30 CONTINUE
339*
340 IF( rightv ) THEN
341*
342* Compute right eigenvectors.
343*
344 is = m
345 DO 80 ki = n, 1, -1
346*
347 IF( somev ) THEN
348 IF( .NOT.SELECT( ki ) )
349 \$ GO TO 80
350 END IF
351 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
352*
353 work( 1 ) = cmone
354*
355* Form right-hand side.
356*
357 DO 40 k = 1, ki - 1
358 work( k ) = -t( k, ki )
359 40 CONTINUE
360*
361* Solve the triangular system:
362* (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
363*
364 DO 50 k = 1, ki - 1
365 t( k, k ) = t( k, k ) - t( ki, ki )
366 IF( cabs1( t( k, k ) ).LT.smin )
367 \$ t( k, k ) = smin
368 50 CONTINUE
369*
370 IF( ki.GT.1 ) THEN
371 CALL clatrs( 'Upper', 'No transpose', 'Non-unit', 'Y',
372 \$ ki-1, t, ldt, work( 1 ), scale, rwork,
373 \$ info )
374 work( ki ) = scale
375 END IF
376*
377* Copy the vector x or Q*x to VR and normalize.
378*
379 IF( .NOT.over ) THEN
380 CALL ccopy( ki, work( 1 ), 1, vr( 1, is ), 1 )
381*
382 ii = icamax( ki, vr( 1, is ), 1 )
383 remax = one / cabs1( vr( ii, is ) )
384 CALL csscal( ki, remax, vr( 1, is ), 1 )
385*
386 DO 60 k = ki + 1, n
387 vr( k, is ) = cmzero
388 60 CONTINUE
389 ELSE
390 IF( ki.GT.1 )
391 \$ CALL cgemv( 'N', n, ki-1, cmone, vr, ldvr, work( 1 ),
392 \$ 1, cmplx( scale ), vr( 1, ki ), 1 )
393*
394 ii = icamax( n, vr( 1, ki ), 1 )
395 remax = one / cabs1( vr( ii, ki ) )
396 CALL csscal( n, remax, vr( 1, ki ), 1 )
397 END IF
398*
399* Set back the original diagonal elements of T.
400*
401 DO 70 k = 1, ki - 1
402 t( k, k ) = work( k+n )
403 70 CONTINUE
404*
405 is = is - 1
406 80 CONTINUE
407 END IF
408*
409 IF( leftv ) THEN
410*
411* Compute left eigenvectors.
412*
413 is = 1
414 DO 130 ki = 1, n
415*
416 IF( somev ) THEN
417 IF( .NOT.SELECT( ki ) )
418 \$ GO TO 130
419 END IF
420 smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
421*
422 work( n ) = cmone
423*
424* Form right-hand side.
425*
426 DO 90 k = ki + 1, n
427 work( k ) = -conjg( t( ki, k ) )
428 90 CONTINUE
429*
430* Solve the triangular system:
431* (T(KI+1:N,KI+1:N) - T(KI,KI))**H*X = SCALE*WORK.
432*
433 DO 100 k = ki + 1, n
434 t( k, k ) = t( k, k ) - t( ki, ki )
435 IF( cabs1( t( k, k ) ).LT.smin )
436 \$ t( k, k ) = smin
437 100 CONTINUE
438*
439 IF( ki.LT.n ) THEN
440 CALL clatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
441 \$ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
442 \$ work( ki+1 ), scale, rwork, info )
443 work( ki ) = scale
444 END IF
445*
446* Copy the vector x or Q*x to VL and normalize.
447*
448 IF( .NOT.over ) THEN
449 CALL ccopy( n-ki+1, work( ki ), 1, vl( ki, is ), 1 )
450*
451 ii = icamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
452 remax = one / cabs1( vl( ii, is ) )
453 CALL csscal( n-ki+1, remax, vl( ki, is ), 1 )
454*
455 DO 110 k = 1, ki - 1
456 vl( k, is ) = cmzero
457 110 CONTINUE
458 ELSE
459 IF( ki.LT.n )
460 \$ CALL cgemv( 'N', n, n-ki, cmone, vl( 1, ki+1 ), ldvl,
461 \$ work( ki+1 ), 1, cmplx( scale ),
462 \$ vl( 1, ki ), 1 )
463*
464 ii = icamax( n, vl( 1, ki ), 1 )
465 remax = one / cabs1( vl( ii, ki ) )
466 CALL csscal( n, remax, vl( 1, ki ), 1 )
467 END IF
468*
469* Set back the original diagonal elements of T.
470*
471 DO 120 k = ki + 1, n
472 t( k, k ) = work( k+n )
473 120 CONTINUE
474*
475 is = is + 1
476 130 CONTINUE
477 END IF
478*
479 RETURN
480*
481* End of CTREVC
482*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:71
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine clatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition: clatrs.f:239
real function scasum(N, CX, INCX)
SCASUM
Definition: scasum.f:72
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: