LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
subroutine ctrevc3 ( 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,
integer  LWORK,
real, dimension( * )  RWORK,
integer  LRWORK,
integer  INFO 
)

CTREVC3

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

Purpose:
 CTREVC3 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.

 This uses a Level 3 BLAS version of the back transformation.
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 (MAX(1,LWORK))
[in]LWORK
          LWORK is INTEGER
          The dimension of array WORK. LWORK >= max(1,2*N).
          For optimum performance, LWORK >= N + 2*N*NB, where NB is
          the optimal blocksize.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is REAL array, dimension (LRWORK)
[in]LRWORK
          LRWORK is INTEGER
          The dimension of array RWORK. LRWORK >= max(1,N).

          If LRWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the RWORK array, returns
          this value as the first entry of the RWORK array, and no error
          message related to LRWORK is issued by XERBLA.
[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.
Date
December 2016
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 248 of file ctrevc3.f.

248  IMPLICIT NONE
249 *
250 * -- LAPACK computational routine (version 3.7.0) --
251 * -- LAPACK is a software package provided by Univ. of Tennessee, --
252 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253 * December 2016
254 *
255 * .. Scalar Arguments ..
256  CHARACTER howmny, side
257  INTEGER info, ldt, ldvl, ldvr, lwork, lrwork, m, mm, n
258 * ..
259 * .. Array Arguments ..
260  LOGICAL select( * )
261  REAL rwork( * )
262  COMPLEX t( ldt, * ), vl( ldvl, * ), vr( ldvr, * ),
263  $ work( * )
264 * ..
265 *
266 * =====================================================================
267 *
268 * .. Parameters ..
269  REAL zero, one
270  parameter ( zero = 0.0e+0, one = 1.0e+0 )
271  COMPLEX czero, cone
272  parameter ( czero = ( 0.0e+0, 0.0e+0 ),
273  $ cone = ( 1.0e+0, 0.0e+0 ) )
274  INTEGER nbmin, nbmax
275  parameter ( nbmin = 8, nbmax = 128 )
276 * ..
277 * .. Local Scalars ..
278  LOGICAL allv, bothv, leftv, lquery, over, rightv, somev
279  INTEGER i, ii, is, j, k, ki, iv, maxwrk, nb
280  REAL ovfl, remax, scale, smin, smlnum, ulp, unfl
281  COMPLEX cdum
282 * ..
283 * .. External Functions ..
284  LOGICAL lsame
285  INTEGER ilaenv, icamax
286  REAL slamch, scasum
287  EXTERNAL lsame, ilaenv, icamax, slamch, scasum
288 * ..
289 * .. External Subroutines ..
290  EXTERNAL xerbla, ccopy, claset, csscal, cgemm, cgemv,
291  $ clatrs, slabad
292 * ..
293 * .. Intrinsic Functions ..
294  INTRINSIC abs, REAL, cmplx, conjg, aimag, max
295 * ..
296 * .. Statement Functions ..
297  REAL cabs1
298 * ..
299 * .. Statement Function definitions ..
300  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( cdum ) )
301 * ..
302 * .. Executable Statements ..
303 *
304 * Decode and test the input parameters
305 *
306  bothv = lsame( side, 'B' )
307  rightv = lsame( side, 'R' ) .OR. bothv
308  leftv = lsame( side, 'L' ) .OR. bothv
309 *
310  allv = lsame( howmny, 'A' )
311  over = lsame( howmny, 'B' )
312  somev = lsame( howmny, 'S' )
313 *
314 * Set M to the number of columns required to store the selected
315 * eigenvectors.
316 *
317  IF( somev ) THEN
318  m = 0
319  DO 10 j = 1, n
320  IF( SELECT( j ) )
321  $ m = m + 1
322  10 CONTINUE
323  ELSE
324  m = n
325  END IF
326 *
327  info = 0
328  nb = ilaenv( 1, 'CTREVC', side // howmny, n, -1, -1, -1 )
329  maxwrk = n + 2*n*nb
330  work(1) = maxwrk
331  rwork(1) = n
332  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
333  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
334  info = -1
335  ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
336  info = -2
337  ELSE IF( n.LT.0 ) THEN
338  info = -4
339  ELSE IF( ldt.LT.max( 1, n ) ) THEN
340  info = -6
341  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
342  info = -8
343  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
344  info = -10
345  ELSE IF( mm.LT.m ) THEN
346  info = -11
347  ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
348  info = -14
349  ELSE IF ( lrwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
350  info = -16
351  END IF
352  IF( info.NE.0 ) THEN
353  CALL xerbla( 'CTREVC3', -info )
354  RETURN
355  ELSE IF( lquery ) THEN
356  RETURN
357  END IF
358 *
359 * Quick return if possible.
360 *
361  IF( n.EQ.0 )
362  $ RETURN
363 *
364 * Use blocked version of back-transformation if sufficient workspace.
365 * Zero-out the workspace to avoid potential NaN propagation.
366 *
367  IF( over .AND. lwork .GE. n + 2*n*nbmin ) THEN
368  nb = (lwork - n) / (2*n)
369  nb = min( nb, nbmax )
370  CALL claset( 'F', n, 1+2*nb, czero, czero, work, n )
371  ELSE
372  nb = 1
373  END IF
374 *
375 * Set the constants to control overflow.
376 *
377  unfl = slamch( 'Safe minimum' )
378  ovfl = one / unfl
379  CALL slabad( unfl, ovfl )
380  ulp = slamch( 'Precision' )
381  smlnum = unfl*( n / ulp )
382 *
383 * Store the diagonal elements of T in working array WORK.
384 *
385  DO 20 i = 1, n
386  work( i ) = t( i, i )
387  20 CONTINUE
388 *
389 * Compute 1-norm of each column of strictly upper triangular
390 * part of T to control overflow in triangular solver.
391 *
392  rwork( 1 ) = zero
393  DO 30 j = 2, n
394  rwork( j ) = scasum( j-1, t( 1, j ), 1 )
395  30 CONTINUE
396 *
397  IF( rightv ) THEN
398 *
399 * ============================================================
400 * Compute right eigenvectors.
401 *
402 * IV is index of column in current block.
403 * Non-blocked version always uses IV=NB=1;
404 * blocked version starts with IV=NB, goes down to 1.
405 * (Note the "0-th" column is used to store the original diagonal.)
406  iv = nb
407  is = m
408  DO 80 ki = n, 1, -1
409  IF( somev ) THEN
410  IF( .NOT.SELECT( ki ) )
411  $ GO TO 80
412  END IF
413  smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
414 *
415 * --------------------------------------------------------
416 * Complex right eigenvector
417 *
418  work( ki + iv*n ) = cone
419 *
420 * Form right-hand side.
421 *
422  DO 40 k = 1, ki - 1
423  work( k + iv*n ) = -t( k, ki )
424  40 CONTINUE
425 *
426 * Solve upper triangular system:
427 * [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
428 *
429  DO 50 k = 1, ki - 1
430  t( k, k ) = t( k, k ) - t( ki, ki )
431  IF( cabs1( t( k, k ) ).LT.smin )
432  $ t( k, k ) = smin
433  50 CONTINUE
434 *
435  IF( ki.GT.1 ) THEN
436  CALL clatrs( 'Upper', 'No transpose', 'Non-unit', 'Y',
437  $ ki-1, t, ldt, work( 1 + iv*n ), scale,
438  $ rwork, info )
439  work( ki + iv*n ) = scale
440  END IF
441 *
442 * Copy the vector x or Q*x to VR and normalize.
443 *
444  IF( .NOT.over ) THEN
445 * ------------------------------
446 * no back-transform: copy x to VR and normalize.
447  CALL ccopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
448 *
449  ii = icamax( ki, vr( 1, is ), 1 )
450  remax = one / cabs1( vr( ii, is ) )
451  CALL csscal( ki, remax, vr( 1, is ), 1 )
452 *
453  DO 60 k = ki + 1, n
454  vr( k, is ) = czero
455  60 CONTINUE
456 *
457  ELSE IF( nb.EQ.1 ) THEN
458 * ------------------------------
459 * version 1: back-transform each vector with GEMV, Q*x.
460  IF( ki.GT.1 )
461  $ CALL cgemv( 'N', n, ki-1, cone, vr, ldvr,
462  $ work( 1 + iv*n ), 1, cmplx( scale ),
463  $ vr( 1, ki ), 1 )
464 *
465  ii = icamax( n, vr( 1, ki ), 1 )
466  remax = one / cabs1( vr( ii, ki ) )
467  CALL csscal( n, remax, vr( 1, ki ), 1 )
468 *
469  ELSE
470 * ------------------------------
471 * version 2: back-transform block of vectors with GEMM
472 * zero out below vector
473  DO k = ki + 1, n
474  work( k + iv*n ) = czero
475  END DO
476 *
477 * Columns IV:NB of work are valid vectors.
478 * When the number of vectors stored reaches NB,
479 * or if this was last vector, do the GEMM
480  IF( (iv.EQ.1) .OR. (ki.EQ.1) ) THEN
481  CALL cgemm( 'N', 'N', n, nb-iv+1, ki+nb-iv, cone,
482  $ vr, ldvr,
483  $ work( 1 + (iv)*n ), n,
484  $ czero,
485  $ work( 1 + (nb+iv)*n ), n )
486 * normalize vectors
487  DO k = iv, nb
488  ii = icamax( n, work( 1 + (nb+k)*n ), 1 )
489  remax = one / cabs1( work( ii + (nb+k)*n ) )
490  CALL csscal( n, remax, work( 1 + (nb+k)*n ), 1 )
491  END DO
492  CALL clacpy( 'F', n, nb-iv+1,
493  $ work( 1 + (nb+iv)*n ), n,
494  $ vr( 1, ki ), ldvr )
495  iv = nb
496  ELSE
497  iv = iv - 1
498  END IF
499  END IF
500 *
501 * Restore the original diagonal elements of T.
502 *
503  DO 70 k = 1, ki - 1
504  t( k, k ) = work( k )
505  70 CONTINUE
506 *
507  is = is - 1
508  80 CONTINUE
509  END IF
510 *
511  IF( leftv ) THEN
512 *
513 * ============================================================
514 * Compute left eigenvectors.
515 *
516 * IV is index of column in current block.
517 * Non-blocked version always uses IV=1;
518 * blocked version starts with IV=1, goes up to NB.
519 * (Note the "0-th" column is used to store the original diagonal.)
520  iv = 1
521  is = 1
522  DO 130 ki = 1, n
523 *
524  IF( somev ) THEN
525  IF( .NOT.SELECT( ki ) )
526  $ GO TO 130
527  END IF
528  smin = max( ulp*( cabs1( t( ki, ki ) ) ), smlnum )
529 *
530 * --------------------------------------------------------
531 * Complex left eigenvector
532 *
533  work( ki + iv*n ) = cone
534 *
535 * Form right-hand side.
536 *
537  DO 90 k = ki + 1, n
538  work( k + iv*n ) = -conjg( t( ki, k ) )
539  90 CONTINUE
540 *
541 * Solve conjugate-transposed triangular system:
542 * [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
543 *
544  DO 100 k = ki + 1, n
545  t( k, k ) = t( k, k ) - t( ki, ki )
546  IF( cabs1( t( k, k ) ).LT.smin )
547  $ t( k, k ) = smin
548  100 CONTINUE
549 *
550  IF( ki.LT.n ) THEN
551  CALL clatrs( 'Upper', 'Conjugate transpose', 'Non-unit',
552  $ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
553  $ work( ki+1 + iv*n ), scale, rwork, info )
554  work( ki + iv*n ) = scale
555  END IF
556 *
557 * Copy the vector x or Q*x to VL and normalize.
558 *
559  IF( .NOT.over ) THEN
560 * ------------------------------
561 * no back-transform: copy x to VL and normalize.
562  CALL ccopy( n-ki+1, work( ki + iv*n ), 1, vl(ki,is), 1 )
563 *
564  ii = icamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
565  remax = one / cabs1( vl( ii, is ) )
566  CALL csscal( n-ki+1, remax, vl( ki, is ), 1 )
567 *
568  DO 110 k = 1, ki - 1
569  vl( k, is ) = czero
570  110 CONTINUE
571 *
572  ELSE IF( nb.EQ.1 ) THEN
573 * ------------------------------
574 * version 1: back-transform each vector with GEMV, Q*x.
575  IF( ki.LT.n )
576  $ CALL cgemv( 'N', n, n-ki, cone, vl( 1, ki+1 ), ldvl,
577  $ work( ki+1 + iv*n ), 1, cmplx( scale ),
578  $ vl( 1, ki ), 1 )
579 *
580  ii = icamax( n, vl( 1, ki ), 1 )
581  remax = one / cabs1( vl( ii, ki ) )
582  CALL csscal( n, remax, vl( 1, ki ), 1 )
583 *
584  ELSE
585 * ------------------------------
586 * version 2: back-transform block of vectors with GEMM
587 * zero out above vector
588 * could go from KI-NV+1 to KI-1
589  DO k = 1, ki - 1
590  work( k + iv*n ) = czero
591  END DO
592 *
593 * Columns 1:IV of work are valid vectors.
594 * When the number of vectors stored reaches NB,
595 * or if this was last vector, do the GEMM
596  IF( (iv.EQ.nb) .OR. (ki.EQ.n) ) THEN
597  CALL cgemm( 'N', 'N', n, iv, n-ki+iv, cone,
598  $ vl( 1, ki-iv+1 ), ldvl,
599  $ work( ki-iv+1 + (1)*n ), n,
600  $ czero,
601  $ work( 1 + (nb+1)*n ), n )
602 * normalize vectors
603  DO k = 1, iv
604  ii = icamax( n, work( 1 + (nb+k)*n ), 1 )
605  remax = one / cabs1( work( ii + (nb+k)*n ) )
606  CALL csscal( n, remax, work( 1 + (nb+k)*n ), 1 )
607  END DO
608  CALL clacpy( 'F', n, iv,
609  $ work( 1 + (nb+1)*n ), n,
610  $ vl( 1, ki-iv+1 ), ldvl )
611  iv = 1
612  ELSE
613  iv = iv + 1
614  END IF
615  END IF
616 *
617 * Restore the original diagonal elements of T.
618 *
619  DO 120 k = ki + 1, n
620  t( k, k ) = work( k )
621  120 CONTINUE
622 *
623  is = is + 1
624  130 CONTINUE
625  END IF
626 *
627  RETURN
628 *
629 * End of CTREVC3
630 *
integer function icamax(N, CX, INCX)
ICAMAX
Definition: icamax.f:53
real function scasum(N, CX, INCX)
SCASUM
Definition: scasum.f:54
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
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:241
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: