 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ 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: