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

◆ ztrevc3()

subroutine ztrevc3 ( character  side,
character  howmny,
logical, dimension( * )  select,
integer  n,
complex*16, dimension( ldt, * )  t,
integer  ldt,
complex*16, dimension( ldvl, * )  vl,
integer  ldvl,
complex*16, dimension( ldvr, * )  vr,
integer  ldvr,
integer  mm,
integer  m,
complex*16, dimension( * )  work,
integer  lwork,
double precision, dimension( * )  rwork,
integer  lrwork,
integer  info 
)

ZTREVC3

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

Purpose:
 ZTREVC3 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 ZHSEQR.

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