LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ctrevc3()

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.
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 242 of file ctrevc3.f.

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