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

Download CTREVC + dependencies [TGZ] [ZIP] [TXT]

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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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 slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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: