LAPACK 3.12.1
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 240 of file ztrevc3.f.

243 IMPLICIT NONE
244*
245* -- LAPACK computational routine --
246* -- LAPACK is a software package provided by Univ. of Tennessee, --
247* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
248*
249* .. Scalar Arguments ..
250 CHARACTER HOWMNY, SIDE
251 INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
252* ..
253* .. Array Arguments ..
254 LOGICAL SELECT( * )
255 DOUBLE PRECISION RWORK( * )
256 COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
257 $ WORK( * )
258* ..
259*
260* =====================================================================
261*
262* .. Parameters ..
263 DOUBLE PRECISION ZERO, ONE
264 parameter( zero = 0.0d+0, one = 1.0d+0 )
265 COMPLEX*16 CZERO, CONE
266 parameter( czero = ( 0.0d+0, 0.0d+0 ),
267 $ cone = ( 1.0d+0, 0.0d+0 ) )
268 INTEGER NBMIN, NBMAX
269 parameter( nbmin = 8, nbmax = 128 )
270* ..
271* .. Local Scalars ..
272 LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
273 INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
274 DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
275 COMPLEX*16 CDUM
276* ..
277* .. External Functions ..
278 LOGICAL LSAME
279 INTEGER ILAENV, IZAMAX
280 DOUBLE PRECISION DLAMCH, DZASUM
281 EXTERNAL lsame, ilaenv, izamax, dlamch,
282 $ dzasum
283* ..
284* .. External Subroutines ..
285 EXTERNAL xerbla, zcopy, zdscal, zgemv,
286 $ zlatrs,
288* ..
289* .. Intrinsic Functions ..
290 INTRINSIC abs, dble, dcmplx, conjg, dimag, max
291* ..
292* .. Statement Functions ..
293 DOUBLE PRECISION CABS1
294* ..
295* .. Statement Function definitions ..
296 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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, 'ZTREVC', side // howmny, n, -1, -1, -1 )
325 maxwrk = max( 1, n + 2*n*nb )
326 work(1) = 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( 'ZTREVC3', -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 zlaset( '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 = dlamch( 'Safe minimum' )
374 ovfl = one / unfl
375 ulp = dlamch( '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 ) = dzasum( 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 zlatrs( '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 zcopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
443*
444 ii = izamax( ki, vr( 1, is ), 1 )
445 remax = one / cabs1( vr( ii, is ) )
446 CALL zdscal( 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 zgemv( 'N', n, ki-1, cone, vr, ldvr,
457 $ work( 1 + iv*n ), 1, dcmplx( scale ),
458 $ vr( 1, ki ), 1 )
459*
460 ii = izamax( n, vr( 1, ki ), 1 )
461 remax = one / cabs1( vr( ii, ki ) )
462 CALL zdscal( 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 zgemm( '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 = izamax( n, work( 1 + (nb+k)*n ), 1 )
484 remax = one / cabs1( work( ii + (nb+k)*n ) )
485 CALL zdscal( n, remax, work( 1 + (nb+k)*n ), 1 )
486 END DO
487 CALL zlacpy( '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 zlatrs( 'Upper', 'Conjugate transpose',
547 $ 'Non-unit',
548 $ 'Y', n-ki, t( ki+1, ki+1 ), ldt,
549 $ work( ki+1 + iv*n ), scale, rwork, info )
550 work( ki + iv*n ) = scale
551 END IF
552*
553* Copy the vector x or Q*x to VL and normalize.
554*
555 IF( .NOT.over ) THEN
556* ------------------------------
557* no back-transform: copy x to VL and normalize.
558 CALL zcopy( n-ki+1, work( ki + iv*n ), 1, vl(ki,is),
559 $ 1 )
560*
561 ii = izamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
562 remax = one / cabs1( vl( ii, is ) )
563 CALL zdscal( n-ki+1, remax, vl( ki, is ), 1 )
564*
565 DO 110 k = 1, ki - 1
566 vl( k, is ) = czero
567 110 CONTINUE
568*
569 ELSE IF( nb.EQ.1 ) THEN
570* ------------------------------
571* version 1: back-transform each vector with GEMV, Q*x.
572 IF( ki.LT.n )
573 $ CALL zgemv( 'N', n, n-ki, cone, vl( 1, ki+1 ),
574 $ ldvl,
575 $ work( ki+1 + iv*n ), 1, dcmplx( scale ),
576 $ vl( 1, ki ), 1 )
577*
578 ii = izamax( n, vl( 1, ki ), 1 )
579 remax = one / cabs1( vl( ii, ki ) )
580 CALL zdscal( n, remax, vl( 1, ki ), 1 )
581*
582 ELSE
583* ------------------------------
584* version 2: back-transform block of vectors with GEMM
585* zero out above vector
586* could go from KI-NV+1 to KI-1
587 DO k = 1, ki - 1
588 work( k + iv*n ) = czero
589 END DO
590*
591* Columns 1:IV of work are valid vectors.
592* When the number of vectors stored reaches NB,
593* or if this was last vector, do the GEMM
594 IF( (iv.EQ.nb) .OR. (ki.EQ.n) ) THEN
595 CALL zgemm( 'N', 'N', n, iv, n-ki+iv, cone,
596 $ vl( 1, ki-iv+1 ), ldvl,
597 $ work( ki-iv+1 + (1)*n ), n,
598 $ czero,
599 $ work( 1 + (nb+1)*n ), n )
600* normalize vectors
601 DO k = 1, iv
602 ii = izamax( n, work( 1 + (nb+k)*n ), 1 )
603 remax = one / cabs1( work( ii + (nb+k)*n ) )
604 CALL zdscal( n, remax, work( 1 + (nb+k)*n ), 1 )
605 END DO
606 CALL zlacpy( 'F', n, iv,
607 $ work( 1 + (nb+1)*n ), n,
608 $ vl( 1, ki-iv+1 ), ldvl )
609 iv = 1
610 ELSE
611 iv = iv + 1
612 END IF
613 END IF
614*
615* Restore the original diagonal elements of T.
616*
617 DO 120 k = ki + 1, n
618 t( k, k ) = work( k )
619 120 CONTINUE
620*
621 is = is + 1
622 130 CONTINUE
623 END IF
624*
625 RETURN
626*
627* End of ZTREVC3
628*
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:160
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
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:104
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:238
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: