summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJulie <julie@cs.utk.edu>2016-06-12 23:21:04 -0700
committerJulie <julie@cs.utk.edu>2016-06-12 23:21:04 -0700
commited2ea1af894955ddd1ddfd0acb15e1c07d459f1e (patch)
tree98e082131f1ca4ec697e7522ce524b2e59ca3b24
parentf22614a1a00c722ee0c570a4e3d36af4f1cb2cb6 (diff)
downloadlapack-ed2ea1af894955ddd1ddfd0acb15e1c07d459f1e.tar.gz
lapack-ed2ea1af894955ddd1ddfd0acb15e1c07d459f1e.tar.bz2
lapack-ed2ea1af894955ddd1ddfd0acb15e1c07d459f1e.zip
blocked back-transformation for the non-symmetric eigenvalue problem - Contribution from Mark Gates (UTK)
From mark: It blocks NB gemv calls into one gemm call inside trevc. To do that, it needs a new routine, trevc3, because unfortunately the lwork was not passed into trevc. (I highly recommend all new routines always pass lwork and lrwork, where applicable, to enable future upgrades & to catch lwork bugs.)
Diffstat
-rw-r--r--SRC/CMakeLists.txt8
-rw-r--r--SRC/Makefile8
-rw-r--r--SRC/cgeev.f31
-rw-r--r--SRC/cgeevx.f26
-rw-r--r--SRC/ctrevc3.f630
-rw-r--r--SRC/dgeev.f33
-rw-r--r--SRC/dgeevx.f29
-rw-r--r--SRC/dtrevc3.f1303
-rw-r--r--SRC/ilaenv.f6
-rw-r--r--SRC/sgeev.f33
-rw-r--r--SRC/sgeevx.f30
-rw-r--r--SRC/strevc3.f1303
-rw-r--r--SRC/zgeev.f30
-rw-r--r--SRC/zgeevx.f26
-rw-r--r--SRC/ztrevc3.f630
15 files changed, 4045 insertions, 81 deletions
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index 03441b94..4857f474 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -141,7 +141,7 @@ set(SLASRC
stbrfs.f stbtrs.f stgevc.f stgex2.f stgexc.f stgsen.f
stgsja.f stgsna.f stgsy2.f stgsyl.f stpcon.f stprfs.f stptri.f
stptrs.f
- strcon.f strevc.f strexc.f strrfs.f strsen.f strsna.f strsyl.f
+ strcon.f strevc.f strevc3.f strexc.f strrfs.f strsen.f strsna.f strsyl.f
strti2.f strtri.f strtrs.f stzrzf.f sstemr.f
slansf.f spftrf.f spftri.f spftrs.f ssfrk.f stfsm.f stftri.f stfttp.f
stfttr.f stpttf.f stpttr.f strttf.f strttp.f
@@ -221,7 +221,7 @@ set(CLASRC
ctbcon.f ctbrfs.f ctbtrs.f ctgevc.f ctgex2.f
ctgexc.f ctgsen.f ctgsja.f ctgsna.f ctgsy2.f ctgsyl.f ctpcon.f
ctprfs.f ctptri.f
- ctptrs.f ctrcon.f ctrevc.f ctrexc.f ctrrfs.f ctrsen.f ctrsna.f
+ ctptrs.f ctrcon.f ctrevc.f ctrevc3.f ctrexc.f ctrrfs.f ctrsen.f ctrsna.f
ctrsyl.f ctrti2.f ctrtri.f ctrtrs.f ctzrzf.f cung2l.f cung2r.f
cungbr.f cunghr.f cungl2.f cunglq.f cungql.f cungqr.f cungr2.f
cungrq.f cungtr.f cunm2l.f cunm2r.f cunmbr.f cunmhr.f cunml2.f cunm22.f
@@ -302,7 +302,7 @@ set(DLASRC
dtbrfs.f dtbtrs.f dtgevc.f dtgex2.f dtgexc.f dtgsen.f
dtgsja.f dtgsna.f dtgsy2.f dtgsyl.f dtpcon.f dtprfs.f dtptri.f
dtptrs.f
- dtrcon.f dtrevc.f dtrexc.f dtrrfs.f dtrsen.f dtrsna.f dtrsyl.f
+ dtrcon.f dtrevc.f dtrevc3.f dtrexc.f dtrrfs.f dtrsen.f dtrsna.f dtrsyl.f
dtrti2.f dtrtri.f dtrtrs.f dtzrzf.f dstemr.f
dsgesv.f dsposv.f dlag2s.f slag2d.f dlat2s.f
dlansf.f dpftrf.f dpftri.f dpftrs.f dsfrk.f dtfsm.f dtftri.f dtfttp.f
@@ -383,7 +383,7 @@ set(ZLASRC
ztbcon.f ztbrfs.f ztbtrs.f ztgevc.f ztgex2.f
ztgexc.f ztgsen.f ztgsja.f ztgsna.f ztgsy2.f ztgsyl.f ztpcon.f
ztprfs.f ztptri.f
- ztptrs.f ztrcon.f ztrevc.f ztrexc.f ztrrfs.f ztrsen.f ztrsna.f
+ ztptrs.f ztrcon.f ztrevc.f ztrevc3.f ztrexc.f ztrrfs.f ztrsen.f ztrsna.f
ztrsyl.f ztrti2.f ztrtri.f ztrtrs.f ztzrzf.f zung2l.f
zung2r.f zungbr.f zunghr.f zungl2.f zunglq.f zungql.f zungqr.f zungr2.f
zungrq.f zungtr.f zunm2l.f zunm2r.f zunmbr.f zunmhr.f zunml2.f zunm22.f
diff --git a/SRC/Makefile b/SRC/Makefile
index f2c67e01..d0163dc1 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -150,7 +150,7 @@ SLASRC = \
stbrfs.o stbtrs.o stgevc.o stgex2.o stgexc.o stgsen.o \
stgsja.o stgsna.o stgsy2.o stgsyl.o stpcon.o stprfs.o stptri.o \
stptrs.o \
- strcon.o strevc.o strexc.o strrfs.o strsen.o strsna.o strsyl.o \
+ strcon.o strevc.o strevc3.o strexc.o strrfs.o strsen.o strsna.o strsyl.o \
strti2.o strtri.o strtrs.o stzrzf.o sstemr.o \
slansf.o spftrf.o spftri.o spftrs.o ssfrk.o stfsm.o stftri.o stfttp.o \
stfttr.o stpttf.o stpttr.o strttf.o strttp.o \
@@ -231,7 +231,7 @@ CLASRC = \
ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \
ctgexc.o ctgsen.o ctgsja.o ctgsna.o ctgsy2.o ctgsyl.o ctpcon.o \
ctprfs.o ctptri.o \
- ctptrs.o ctrcon.o ctrevc.o ctrexc.o ctrrfs.o ctrsen.o ctrsna.o \
+ ctptrs.o ctrcon.o ctrevc.o ctrevc3.o ctrexc.o ctrrfs.o ctrsen.o ctrsna.o \
ctrsyl.o ctrti2.o ctrtri.o ctrtrs.o ctzrzf.o cung2l.o cung2r.o \
cungbr.o cunghr.o cungl2.o cunglq.o cungql.o cungqr.o cungr2.o \
cungrq.o cungtr.o cunm2l.o cunm2r.o cunmbr.o cunmhr.o cunml2.o cunm22.o \
@@ -316,7 +316,7 @@ DLASRC = \
dtbcon.o dtbrfs.o dtbtrs.o dtgevc.o dtgex2.o dtgexc.o dtgsen.o \
dtgsja.o dtgsna.o dtgsy2.o dtgsyl.o dtpcon.o dtprfs.o dtptri.o \
dtptrs.o \
- dtrcon.o dtrevc.o dtrexc.o dtrrfs.o dtrsen.o dtrsna.o dtrsyl.o \
+ dtrcon.o dtrevc.o dtrevc3.o dtrexc.o dtrrfs.o dtrsen.o dtrsna.o dtrsyl.o \
dtrti2.o dtrtri.o dtrtrs.o dtzrzf.o dstemr.o \
dsgesv.o dsposv.o dlag2s.o slag2d.o dlat2s.o \
dlansf.o dpftrf.o dpftri.o dpftrs.o dsfrk.o dtfsm.o dtftri.o dtfttp.o \
@@ -400,7 +400,7 @@ ZLASRC = \
ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \
ztgexc.o ztgsen.o ztgsja.o ztgsna.o ztgsy2.o ztgsyl.o ztpcon.o \
ztprfs.o ztptri.o \
- ztptrs.o ztrcon.o ztrevc.o ztrexc.o ztrrfs.o ztrsen.o ztrsna.o \
+ ztptrs.o ztrcon.o ztrevc.o ztrevc3.o ztrexc.o ztrrfs.o ztrsen.o ztrsna.o \
ztrsyl.o ztrti2.o ztrtri.o ztrtrs.o ztzrzf.o zung2l.o \
zung2r.o zungbr.o zunghr.o zungl2.o zunglq.o zungql.o zungqr.o zungr2.o \
zungrq.o zungtr.o zunm2l.o zunm2r.o zunmbr.o zunmhr.o zunml2.o zunm22.o \
diff --git a/SRC/cgeev.f b/SRC/cgeev.f
index 0f48322a..a888c64f 100644
--- a/SRC/cgeev.f
+++ b/SRC/cgeev.f
@@ -176,6 +176,7 @@
* =====================================================================
SUBROUTINE CGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
$ WORK, LWORK, RWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -202,7 +203,7 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
CHARACTER SIDE
INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
- $ IWRK, K, MAXWRK, MINWRK, NOUT
+ $ IWRK, K, LWORK_TREVC, MAXWRK, MINWRK, NOUT
REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
COMPLEX TMP
* ..
@@ -212,7 +213,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
- $ CSCAL, CSSCAL, CTREVC, CUNGHR, SLABAD, XERBLA
+ $ CSCAL, CSSCAL, CTREVC3, CUNGHR, SLABAD, XERBLA
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -244,7 +245,6 @@
ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
INFO = -10
END IF
-
*
* Compute workspace
* (Note: Comments in the code beginning "Workspace:" describe the
@@ -267,18 +267,28 @@
IF( WANTVL ) THEN
MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
$ ' ', N, 1, N, -1 ) )
+ CALL CTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
$ ' ', N, 1, N, -1 ) )
+ CALL CTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
ELSE
CALL CHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
END IF
WORK( 1 ) = MAXWRK
@@ -413,12 +423,13 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (CWorkspace: need 2*N)
+* (CWorkspace: need 2*N, prefer N + 2*N*NB)
* (RWorkspace: need 2*N)
*
IRWORK = IBAL + N
- CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
+ CALL CTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
+ $ RWORK( IRWORK ), N, IERR )
END IF
*
IF( WANTVL ) THEN
diff --git a/SRC/cgeevx.f b/SRC/cgeevx.f
index 1a80e8c4..b62f070c 100644
--- a/SRC/cgeevx.f
+++ b/SRC/cgeevx.f
@@ -284,6 +284,7 @@
SUBROUTINE CGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
$ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
$ RCONDV, WORK, LWORK, RWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -312,8 +313,8 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
$ WNTSNN, WNTSNV
CHARACTER JOB, SIDE
- INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
- $ MINWRK, NOUT
+ INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
REAL ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
COMPLEX TMP
* ..
@@ -323,7 +324,7 @@
* ..
* .. External Subroutines ..
EXTERNAL CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
- $ CSCAL, CSSCAL, CTREVC, CTRSNA, CUNGHR, SLABAD,
+ $ CSCAL, CSSCAL, CTREVC3, CTRSNA, CUNGHR, SLABAD,
$ SLASCL, XERBLA
* ..
* .. External Functions ..
@@ -387,9 +388,19 @@
MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
*
IF( WANTVL ) THEN
+ CALL CTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, LWORK_TREVC )
CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
$ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
+ CALL CTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, LWORK_TREVC )
CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
$ WORK, -1, INFO )
ELSE
@@ -401,7 +412,7 @@
$ WORK, -1, INFO )
END IF
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
*
IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
MINWRK = 2*N
@@ -567,11 +578,12 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (CWorkspace: need 2*N)
+* (CWorkspace: need 2*N, prefer N + 2*N*NB)
* (RWorkspace: need N)
*
- CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), RWORK, IERR )
+ CALL CTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
+ $ RWORK, N, IERR )
END IF
*
* Compute condition numbers if desired
diff --git a/SRC/ctrevc3.f b/SRC/ctrevc3.f
new file mode 100644
index 00000000..00d3b946
--- /dev/null
+++ b/SRC/ctrevc3.f
@@ -0,0 +1,630 @@
+*> \brief \b CTREVC3
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download CTREVC3 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrevc3.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrevc3.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrevc3.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE CTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+* VR, LDVR, MM, M, WORK, LWORK, RWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER HOWMNY, SIDE
+* INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* REAL RWORK( * )
+* COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+* $ WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*> SIDE is CHARACTER*1
+*> = 'R': compute right eigenvectors only;
+*> = 'L': compute left eigenvectors only;
+*> = 'B': compute both right and left eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] HOWMNY
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in] SELECT
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix T. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] T
+*> \verbatim
+*> T is COMPLEX array, dimension (LDT,N)
+*> The upper triangular matrix T. T is modified, but restored
+*> on exit.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*> LDT is INTEGER
+*> The leading dimension of the array T. LDT >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] VL
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the array VL.
+*> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[in,out] VR
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the array VR.
+*> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[in] MM
+*> \verbatim
+*> MM is INTEGER
+*> The number of columns in the arrays VL and/or VR. MM >= M.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is REAL array, dimension (LRWORK)
+*> \endverbatim
+*>
+*> \param[in] LRWORK
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+* @generated from ztrevc3.f, fortran z -> c, Tue Apr 19 01:47:44 2016
+*
+*> \ingroup complexOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> 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|.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE CTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
+ $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
+ IMPLICIT NONE
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER HOWMNY, SIDE
+ INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+ LOGICAL SELECT( * )
+ REAL RWORK( * )
+ COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+ $ WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ COMPLEX CZERO, CONE
+ PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
+ $ CONE = ( 1.0E+0, 0.0E+0 ) )
+ INTEGER NBMIN, NBMAX
+ PARAMETER ( NBMIN = 8, NBMAX = 128 )
+* ..
+* .. Local Scalars ..
+ LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
+ INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
+ REAL OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
+ COMPLEX CDUM
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ILAENV, ICAMAX
+ REAL SLAMCH, SCASUM
+ EXTERNAL LSAME, ILAENV, ICAMAX, SLAMCH, SCASUM
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, CCOPY, CSSCAL, CGEMV, CLATRS
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, REAL, CMPLX, CONJG, AIMAG, MAX
+* ..
+* .. Statement Functions ..
+ REAL CABS1
+* ..
+* .. Statement Function definitions ..
+ CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
+* ..
+* .. Executable Statements ..
+*
+* Decode and test the input parameters
+*
+ BOTHV = LSAME( SIDE, 'B' )
+ RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
+ LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
+*
+ ALLV = LSAME( HOWMNY, 'A' )
+ OVER = LSAME( HOWMNY, 'B' )
+ SOMEV = LSAME( HOWMNY, 'S' )
+*
+* Set M to the number of columns required to store the selected
+* eigenvectors.
+*
+ IF( SOMEV ) THEN
+ M = 0
+ DO 10 J = 1, N
+ IF( SELECT( J ) )
+ $ M = M + 1
+ 10 CONTINUE
+ ELSE
+ M = N
+ END IF
+*
+ INFO = 0
+ NB = ILAENV( 1, 'CTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
+ MAXWRK = N + 2*N*NB
+ WORK(1) = MAXWRK
+ RWORK(1) = N
+ LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
+ IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
+ INFO = -1
+ ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
+ INFO = -8
+ ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
+ INFO = -10
+ ELSE IF( MM.LT.M ) THEN
+ INFO = -11
+ ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -14
+ ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -16
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'CTREVC3', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Use blocked version of back-transformation if sufficient workspace.
+* Zero-out the workspace to avoid potential NaN propagation.
+*
+ IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
+ NB = (LWORK - N) / (2*N)
+ NB = MIN( NB, NBMAX )
+ CALL CLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
+ ELSE
+ NB = 1
+ END IF
+*
+* Set the constants to control overflow.
+*
+ UNFL = SLAMCH( 'Safe minimum' )
+ OVFL = ONE / UNFL
+ CALL SLABAD( UNFL, OVFL )
+ ULP = SLAMCH( 'Precision' )
+ SMLNUM = UNFL*( N / ULP )
+*
+* Store the diagonal elements of T in working array WORK.
+*
+ DO 20 I = 1, N
+ WORK( I ) = T( I, I )
+ 20 CONTINUE
+*
+* Compute 1-norm of each column of strictly upper triangular
+* part of T to control overflow in triangular solver.
+*
+ RWORK( 1 ) = ZERO
+ DO 30 J = 2, N
+ RWORK( J ) = SCASUM( J-1, T( 1, J ), 1 )
+ 30 CONTINUE
+*
+ IF( RIGHTV ) THEN
+*
+* ============================================================
+* Compute right eigenvectors.
+*
+* IV is index of column in current block.
+* Non-blocked version always uses IV=NB=1;
+* blocked version starts with IV=NB, goes down to 1.
+* (Note the "0-th" column is used to store the original diagonal.)
+ IV = NB
+ IS = M
+ DO 80 KI = N, 1, -1
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 80
+ END IF
+ SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
+*
+* --------------------------------------------------------
+* Complex right eigenvector
+*
+ WORK( KI + IV*N ) = CONE
+*
+* Form right-hand side.
+*
+ DO 40 K = 1, KI - 1
+ WORK( K + IV*N ) = -T( K, KI )
+ 40 CONTINUE
+*
+* Solve upper triangular system:
+* [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
+*
+ DO 50 K = 1, KI - 1
+ T( K, K ) = T( K, K ) - T( KI, KI )
+ IF( CABS1( T( K, K ) ).LT.SMIN )
+ $ T( K, K ) = SMIN
+ 50 CONTINUE
+*
+ IF( KI.GT.1 ) THEN
+ CALL CLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
+ $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
+ $ RWORK, INFO )
+ WORK( KI + IV*N ) = SCALE
+ END IF
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL CCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
+*
+ II = ICAMAX( KI, VR( 1, IS ), 1 )
+ REMAX = ONE / CABS1( VR( II, IS ) )
+ CALL CSSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 60 K = KI + 1, N
+ VR( K, IS ) = CZERO
+ 60 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.1 )
+ $ CALL CGEMV( 'N', N, KI-1, CONE, VR, LDVR,
+ $ WORK( 1 + IV*N ), 1, CMPLX( SCALE ),
+ $ VR( 1, KI ), 1 )
+*
+ II = ICAMAX( N, VR( 1, KI ), 1 )
+ REMAX = ONE / CABS1( VR( II, KI ) )
+ CALL CSSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + IV*N ) = CZERO
+ END DO
+*
+* Columns IV:NB of work are valid vectors.
+* When the number of vectors stored reaches NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
+ CALL CGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
+ $ VR, LDVR,
+ $ WORK( 1 + (IV)*N ), N,
+ $ CZERO,
+ $ WORK( 1 + (NB+IV)*N ), N )
+* normalize vectors
+ DO K = IV, NB
+ II = ICAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
+ CALL CSSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL CLACPY( 'F', N, NB-IV+1,
+ $ WORK( 1 + (NB+IV)*N ), N,
+ $ VR( 1, KI ), LDVR )
+ IV = NB
+ ELSE
+ IV = IV - 1
+ END IF
+ END IF
+*
+* Restore the original diagonal elements of T.
+*
+ DO 70 K = 1, KI - 1
+ T( K, K ) = WORK( K )
+ 70 CONTINUE
+*
+ IS = IS - 1
+ 80 CONTINUE
+ END IF
+*
+ IF( LEFTV ) THEN
+*
+* ============================================================
+* Compute left eigenvectors.
+*
+* IV is index of column in current block.
+* Non-blocked version always uses IV=1;
+* blocked version starts with IV=1, goes up to NB.
+* (Note the "0-th" column is used to store the original diagonal.)
+ IV = 1
+ IS = 1
+ DO 130 KI = 1, N
+*
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 130
+ END IF
+ SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
+*
+* --------------------------------------------------------
+* Complex left eigenvector
+*
+ WORK( KI + IV*N ) = CONE
+*
+* Form right-hand side.
+*
+ DO 90 K = KI + 1, N
+ WORK( K + IV*N ) = -CONJG( T( KI, K ) )
+ 90 CONTINUE
+*
+* Solve conjugate-transposed triangular system:
+* [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
+*
+ DO 100 K = KI + 1, N
+ T( K, K ) = T( K, K ) - T( KI, KI )
+ IF( CABS1( T( K, K ) ).LT.SMIN )
+ $ T( K, K ) = SMIN
+ 100 CONTINUE
+*
+ IF( KI.LT.N ) THEN
+ CALL CLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
+ $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
+ $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
+ WORK( KI + IV*N ) = SCALE
+ END IF
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL CCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
+*
+ II = ICAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
+ REMAX = ONE / CABS1( VL( II, IS ) )
+ CALL CSSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+*
+ DO 110 K = 1, KI - 1
+ VL( K, IS ) = CZERO
+ 110 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N )
+ $ CALL CGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
+ $ WORK( KI+1 + IV*N ), 1, CMPLX( SCALE ),
+ $ VL( 1, KI ), 1 )
+*
+ II = ICAMAX( N, VL( 1, KI ), 1 )
+ REMAX = ONE / CABS1( VL( II, KI ) )
+ CALL CSSCAL( N, REMAX, VL( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + IV*N ) = CZERO
+ END DO
+*
+* Columns 1:IV of work are valid vectors.
+* When the number of vectors stored reaches NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
+ CALL CGEMM( 'N', 'N', N, IV, N-KI+IV, ONE,
+ $ VL( 1, KI-IV+1 ), LDVL,
+ $ WORK( KI-IV+1 + (1)*N ), N,
+ $ CZERO,
+ $ WORK( 1 + (NB+1)*N ), N )
+* normalize vectors
+ DO K = 1, IV
+ II = ICAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
+ CALL CSSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL CLACPY( 'F', N, IV,
+ $ WORK( 1 + (NB+1)*N ), N,
+ $ VL( 1, KI-IV+1 ), LDVL )
+ IV = 1
+ ELSE
+ IV = IV + 1
+ END IF
+ END IF
+*
+* Restore the original diagonal elements of T.
+*
+ DO 120 K = KI + 1, N
+ T( K, K ) = WORK( K )
+ 120 CONTINUE
+*
+ IS = IS + 1
+ 130 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of CTREVC3
+*
+ END
diff --git a/SRC/dgeev.f b/SRC/dgeev.f
index 328eaa39..1c92b7e3 100644
--- a/SRC/dgeev.f
+++ b/SRC/dgeev.f
@@ -188,6 +188,7 @@
* =====================================================================
SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
$ LDVR, WORK, LWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -213,7 +214,7 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
CHARACTER SIDE
INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
- $ MAXWRK, MINWRK, NOUT
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
$ SN
* ..
@@ -223,7 +224,7 @@
* ..
* .. External Subroutines ..
EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
- $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
+ $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
$ XERBLA
* ..
* .. External Functions ..
@@ -279,24 +280,34 @@
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
$ 'DORGHR', ' ', N, 1, N, -1 ) )
CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
+ CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR, N, NOUT,
+ $ WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE IF( WANTVR ) THEN
MINWRK = 4*N
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
$ 'DORGHR', ' ', N, 1, N, -1 ) )
CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
+ CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR, N, NOUT,
+ $ WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE
MINWRK = 3*N
CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
END IF
MAXWRK = MAX( MAXWRK, MINWRK )
@@ -426,10 +437,10 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (Workspace: need 4*N)
+* (Workspace: need 4*N, prefer N + N + 2*N*NB)
*
- CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), IERR )
+ CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
END IF
*
IF( WANTVL ) THEN
diff --git a/SRC/dgeevx.f b/SRC/dgeevx.f
index 8d80d782..d2ba08f0 100644
--- a/SRC/dgeevx.f
+++ b/SRC/dgeevx.f
@@ -302,6 +302,7 @@
SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
$ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
$ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -330,8 +331,8 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
$ WNTSNN, WNTSNV
CHARACTER JOB, SIDE
- INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
- $ MINWRK, NOUT
+ INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
$ SN
* ..
@@ -341,7 +342,7 @@
* ..
* .. External Subroutines ..
EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
- $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
+ $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC3,
$ DTRSNA, XERBLA
* ..
* .. External Functions ..
@@ -366,8 +367,8 @@
WNTSNE = LSAME( SENSE, 'E' )
WNTSNV = LSAME( SENSE, 'V' )
WNTSNB = LSAME( SENSE, 'B' )
- IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC,
- $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
+ IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' )
+ $ .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
$ THEN
INFO = -1
ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
@@ -406,9 +407,19 @@
MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
*
IF( WANTVL ) THEN
+ CALL DTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
$ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
+ CALL DTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
$ WORK, -1, INFO )
ELSE
@@ -420,7 +431,7 @@
$ LDVR, WORK, -1, INFO )
END IF
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
*
IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
MINWRK = 2*N
@@ -580,10 +591,10 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (Workspace: need 3*N)
+* (Workspace: need 3*N, prefer N + 2*N*NB)
*
- CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), IERR )
+ CALL DTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
END IF
*
* Compute condition numbers if desired
diff --git a/SRC/dtrevc3.f b/SRC/dtrevc3.f
new file mode 100644
index 00000000..ba5abb5f
--- /dev/null
+++ b/SRC/dtrevc3.f
@@ -0,0 +1,1303 @@
+*> \brief \b DTREVC3
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DTREVC3 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+* VR, LDVR, MM, M, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER HOWMNY, SIDE
+* INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+* $ WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DTREVC3 computes some or all of the right and/or left eigenvectors of
+*> a real upper quasi-triangular matrix T.
+*> Matrices of this type are produced by the Schur factorization of
+*> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
+*>
+*> The right eigenvector x and the left eigenvector y of T corresponding
+*> to an eigenvalue w are defined by:
+*>
+*> T*x = w*x, (y**T)*T = w*(y**T)
+*>
+*> where y**T denotes the transpose of the vector y.
+*> The eigenvalues are not input to this routine, but are read directly
+*> from the diagonal blocks 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 orthogonal 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*> SIDE is CHARACTER*1
+*> = 'R': compute right eigenvectors only;
+*> = 'L': compute left eigenvectors only;
+*> = 'B': compute both right and left eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] HOWMNY
+*> \verbatim
+*> HOWMNY is CHARACTER*1
+*> = 'A': compute all right and/or left eigenvectors;
+*> = 'B': compute all right and/or left eigenvectors,
+*> backtransformed by the matrices in VR and/or VL;
+*> = 'S': compute selected right and/or left eigenvectors,
+*> as indicated by the logical array SELECT.
+*> \endverbatim
+*>
+*> \param[in,out] SELECT
+*> \verbatim
+*> SELECT is LOGICAL array, dimension (N)
+*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
+*> computed.
+*> If w(j) is a real eigenvalue, the corresponding real
+*> eigenvector is computed if SELECT(j) is .TRUE..
+*> If w(j) and w(j+1) are the real and imaginary parts of a
+*> complex eigenvalue, the corresponding complex eigenvector is
+*> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
+*> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
+*> .FALSE..
+*> Not referenced if HOWMNY = 'A' or 'B'.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix T. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*> T is DOUBLE PRECISION array, dimension (LDT,N)
+*> The upper quasi-triangular matrix T in Schur canonical form.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*> LDT is INTEGER
+*> The leading dimension of the array T. LDT >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] VL
+*> \verbatim
+*> VL is DOUBLE PRECISION 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 orthogonal matrix Q
+*> of Schur vectors returned by DHSEQR).
+*> 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.
+*> A complex eigenvector corresponding to a complex eigenvalue
+*> is stored in two consecutive columns, the first holding the
+*> real part, and the second the imaginary part.
+*> Not referenced if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the array VL.
+*> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[in,out] VR
+*> \verbatim
+*> VR is DOUBLE PRECISION 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 orthogonal matrix Q
+*> of Schur vectors returned by DHSEQR).
+*> 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.
+*> A complex eigenvector corresponding to a complex eigenvalue
+*> is stored in two consecutive columns, the first holding the
+*> real part and the second the imaginary part.
+*> Not referenced if SIDE = 'L'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the array VR.
+*> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[in] MM
+*> \verbatim
+*> MM is INTEGER
+*> The number of columns in the arrays VL and/or VR. MM >= M.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> 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 real eigenvector occupies one column and each
+*> selected complex eigenvector occupies two columns.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of array WORK. LWORK >= max(1,3*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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+* @precisions fortran d -> s
+*
+*> \ingroup doubleOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> 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|.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+ $ VR, LDVR, MM, M, WORK, LWORK, INFO )
+ IMPLICIT NONE
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER HOWMNY, SIDE
+ INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+ LOGICAL SELECT( * )
+ DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+ $ WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ INTEGER NBMIN, NBMAX
+ PARAMETER ( NBMIN = 8, NBMAX = 128 )
+* ..
+* .. Local Scalars ..
+ LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
+ $ RIGHTV, SOMEV
+ INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
+ $ IV, MAXWRK, NB, KI2
+ DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
+ $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
+ $ XNORM
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX, ILAENV
+ DOUBLE PRECISION DDOT, DLAMCH
+ EXTERNAL LSAME, IDAMAX, ILAENV, DDOT, DLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, SQRT
+* ..
+* .. Local Arrays ..
+ DOUBLE PRECISION X( 2, 2 )
+ INTEGER ISCOMPLEX( NBMAX )
+* ..
+* .. Executable Statements ..
+*
+* Decode and test the input parameters
+*
+ BOTHV = LSAME( SIDE, 'B' )
+ RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
+ LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
+*
+ ALLV = LSAME( HOWMNY, 'A' )
+ OVER = LSAME( HOWMNY, 'B' )
+ SOMEV = LSAME( HOWMNY, 'S' )
+*
+ INFO = 0
+ NB = ILAENV( 1, 'DTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
+ MAXWRK = N + 2*N*NB
+ WORK(1) = MAXWRK
+ LQUERY = ( LWORK.EQ.-1 )
+ IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
+ INFO = -1
+ ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
+ INFO = -8
+ ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
+ INFO = -10
+ ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -14
+ ELSE
+*
+* Set M to the number of columns required to store the selected
+* eigenvectors, standardize the array SELECT if necessary, and
+* test MM.
+*
+ IF( SOMEV ) THEN
+ M = 0
+ PAIR = .FALSE.
+ DO 10 J = 1, N
+ IF( PAIR ) THEN
+ PAIR = .FALSE.
+ SELECT( J ) = .FALSE.
+ ELSE
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).EQ.ZERO ) THEN
+ IF( SELECT( J ) )
+ $ M = M + 1
+ ELSE
+ PAIR = .TRUE.
+ IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
+ SELECT( J ) = .TRUE.
+ M = M + 2
+ END IF
+ END IF
+ ELSE
+ IF( SELECT( N ) )
+ $ M = M + 1
+ END IF
+ END IF
+ 10 CONTINUE
+ ELSE
+ M = N
+ END IF
+*
+ IF( MM.LT.M ) THEN
+ INFO = -11
+ END IF
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DTREVC3', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Use blocked version of back-transformation if sufficient workspace.
+* Zero-out the workspace to avoid potential NaN propagation.
+*
+ IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
+ NB = (LWORK - N) / (2*N)
+ NB = MIN( NB, NBMAX )
+ CALL DLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
+ ELSE
+ NB = 1
+ END IF
+*
+* Set the constants to control overflow.
+*
+ UNFL = DLAMCH( 'Safe minimum' )
+ OVFL = ONE / UNFL
+ CALL DLABAD( UNFL, OVFL )
+ ULP = DLAMCH( 'Precision' )
+ SMLNUM = UNFL*( N / ULP )
+ BIGNUM = ( ONE-ULP ) / SMLNUM
+*
+* Compute 1-norm of each column of strictly upper triangular
+* part of T to control overflow in triangular solver.
+*
+ WORK( 1 ) = ZERO
+ DO 30 J = 2, N
+ WORK( J ) = ZERO
+ DO 20 I = 1, J - 1
+ WORK( J ) = WORK( J ) + ABS( T( I, J ) )
+ 20 CONTINUE
+ 30 CONTINUE
+*
+* Index IP is used to specify the real or complex eigenvalue:
+* IP = 0, real eigenvalue,
+* 1, first of conjugate complex pair: (wr,wi)
+* -1, second of conjugate complex pair: (wr,wi)
+* ISCOMPLEX array stores IP for each column in current block.
+*
+ IF( RIGHTV ) THEN
+*
+* ============================================================
+* Compute right eigenvectors.
+*
+* IV is index of column in current block.
+* For complex right vector, uses IV-1 for real part and IV for complex part.
+* Non-blocked version always uses IV=2;
+* blocked version starts with IV=NB, goes down to 1 or 2.
+* (Note the "0-th" column is used for 1-norms computed above.)
+ IV = 2
+ IF( NB.GT.2 ) THEN
+ IV = NB
+ END IF
+
+ IP = 0
+ IS = M
+ DO 140 KI = N, 1, -1
+ IF( IP.EQ.-1 ) THEN
+* previous iteration (ki+1) was second of conjugate pair,
+* so this ki is first of conjugate pair; skip to end of loop
+ IP = 1
+ GO TO 140
+ ELSE IF( KI.EQ.1 ) THEN
+* last column, so this ki must be real eigenvalue
+ IP = 0
+ ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
+* zero on sub-diagonal, so this ki is real eigenvalue
+ IP = 0
+ ELSE
+* non-zero on sub-diagonal, so this ki is second of conjugate pair
+ IP = -1
+ END IF
+
+ IF( SOMEV ) THEN
+ IF( IP.EQ.0 ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 140
+ ELSE
+ IF( .NOT.SELECT( KI-1 ) )
+ $ GO TO 140
+ END IF
+ END IF
+*
+* Compute the KI-th eigenvalue (WR,WI).
+*
+ WR = T( KI, KI )
+ WI = ZERO
+ IF( IP.NE.0 )
+ $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
+ $ SQRT( ABS( T( KI-1, KI ) ) )
+ SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
+*
+ IF( IP.EQ.0 ) THEN
+*
+* --------------------------------------------------------
+* Real right eigenvector
+*
+ WORK( KI + IV*N ) = ONE
+*
+* Form right-hand side.
+*
+ DO 50 K = 1, KI - 1
+ WORK( K + IV*N ) = -T( K, KI )
+ 50 CONTINUE
+*
+* Solve upper quasi-triangular system:
+* [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
+*
+ JNXT = KI - 1
+ DO 60 J = KI - 1, 1, -1
+ IF( J.GT.JNXT )
+ $ GO TO 60
+ J1 = J
+ J2 = J
+ JNXT = J - 1
+ IF( J.GT.1 ) THEN
+ IF( T( J, J-1 ).NE.ZERO ) THEN
+ J1 = J - 1
+ JNXT = J - 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+ CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale X(1,1) to avoid overflow when updating
+* the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
+ WORK( J+IV*N ) = X( 1, 1 )
+*
+* Update right-hand side
+*
+ CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
+ $ WORK( 1+IV*N ), 1 )
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+ CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
+ $ T( J-1, J-1 ), LDT, ONE, ONE,
+ $ WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
+ $ SCALE, XNORM, IERR )
+*
+* Scale X(1,1) and X(2,1) to avoid overflow when
+* updating the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ BETA = MAX( WORK( J-1 ), WORK( J ) )
+ IF( BETA.GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ X( 2, 1 ) = X( 2, 1 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL DSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
+ WORK( J-1+IV*N ) = X( 1, 1 )
+ WORK( J +IV*N ) = X( 2, 1 )
+*
+* Update right-hand side
+*
+ CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
+ $ WORK( 1+IV*N ), 1 )
+ CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
+ $ WORK( 1+IV*N ), 1 )
+ END IF
+ 60 CONTINUE
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL DCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
+*
+ II = IDAMAX( KI, VR( 1, IS ), 1 )
+ REMAX = ONE / ABS( VR( II, IS ) )
+ CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 70 K = KI + 1, N
+ VR( K, IS ) = ZERO
+ 70 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.1 )
+ $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
+ $ WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
+ $ VR( 1, KI ), 1 )
+*
+ II = IDAMAX( N, VR( 1, KI ), 1 )
+ REMAX = ONE / ABS( VR( II, KI ) )
+ CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + IV*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+* back-transform and normalization is done below
+ END IF
+ ELSE
+*
+* --------------------------------------------------------
+* Complex right eigenvector.
+*
+* Initial solve
+* [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
+* [ ( T(KI, KI-1) T(KI, KI) ) ]
+*
+ IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
+ WORK( KI-1 + (IV-1)*N ) = ONE
+ WORK( KI + (IV )*N ) = WI / T( KI-1, KI )
+ ELSE
+ WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
+ WORK( KI + (IV )*N ) = ONE
+ END IF
+ WORK( KI + (IV-1)*N ) = ZERO
+ WORK( KI-1 + (IV )*N ) = ZERO
+*
+* Form right-hand side.
+*
+ DO 80 K = 1, KI - 2
+ WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
+ WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(K,KI )
+ 80 CONTINUE
+*
+* Solve upper quasi-triangular system:
+* [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
+*
+ JNXT = KI - 2
+ DO 90 J = KI - 2, 1, -1
+ IF( J.GT.JNXT )
+ $ GO TO 90
+ J1 = J
+ J2 = J
+ JNXT = J - 1
+ IF( J.GT.1 ) THEN
+ IF( T( J, J-1 ).NE.ZERO ) THEN
+ J1 = J - 1
+ JNXT = J - 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+ CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
+ $ WR, WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale X(1,1) and X(1,2) to avoid overflow when
+* updating the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ X( 1, 2 ) = X( 1, 2 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
+ CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
+ END IF
+ WORK( J+(IV-1)*N ) = X( 1, 1 )
+ WORK( J+(IV )*N ) = X( 1, 2 )
+*
+* Update the right-hand side
+*
+ CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+ CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
+ $ T( J-1, J-1 ), LDT, ONE, ONE,
+ $ WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
+ $ SCALE, XNORM, IERR )
+*
+* Scale X to avoid overflow when updating
+* the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ BETA = MAX( WORK( J-1 ), WORK( J ) )
+ IF( BETA.GT.BIGNUM / XNORM ) THEN
+ REC = ONE / XNORM
+ X( 1, 1 ) = X( 1, 1 )*REC
+ X( 1, 2 ) = X( 1, 2 )*REC
+ X( 2, 1 ) = X( 2, 1 )*REC
+ X( 2, 2 ) = X( 2, 2 )*REC
+ SCALE = SCALE*REC
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL DSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
+ CALL DSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
+ END IF
+ WORK( J-1+(IV-1)*N ) = X( 1, 1 )
+ WORK( J +(IV-1)*N ) = X( 2, 1 )
+ WORK( J-1+(IV )*N ) = X( 1, 2 )
+ WORK( J +(IV )*N ) = X( 2, 2 )
+*
+* Update the right-hand side
+*
+ CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+ CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+ END IF
+ 90 CONTINUE
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL DCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
+ CALL DCOPY( KI, WORK( 1+(IV )*N ), 1, VR(1,IS ), 1 )
+*
+ EMAX = ZERO
+ DO 100 K = 1, KI
+ EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
+ $ ABS( VR( K, IS ) ) )
+ 100 CONTINUE
+ REMAX = ONE / EMAX
+ CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
+ CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 110 K = KI + 1, N
+ VR( K, IS-1 ) = ZERO
+ VR( K, IS ) = ZERO
+ 110 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.2 ) THEN
+ CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
+ $ WORK( 1 + (IV-1)*N ), 1,
+ $ WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
+ CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
+ $ WORK( 1 + (IV)*N ), 1,
+ $ WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
+ ELSE
+ CALL DSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
+ CALL DSCAL( N, WORK(KI +(IV )*N), VR(1,KI ), 1)
+ END IF
+*
+ EMAX = ZERO
+ DO 120 K = 1, N
+ EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
+ $ ABS( VR( K, KI ) ) )
+ 120 CONTINUE
+ REMAX = ONE / EMAX
+ CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
+ CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + (IV-1)*N ) = ZERO
+ WORK( K + (IV )*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV-1 ) = -IP
+ ISCOMPLEX( IV ) = IP
+ IV = IV - 1
+* back-transform and normalization is done below
+ END IF
+ END IF
+
+ IF( NB.GT.1 ) THEN
+* --------------------------------------------------------
+* Blocked version of back-transform
+* For complex case, KI2 includes both vectors (KI-1 and KI)
+ IF( IP.EQ.0 ) THEN
+ KI2 = KI
+ ELSE
+ KI2 = KI - 1
+ END IF
+
+* Columns IV:NB of work are valid vectors.
+* When the number of vectors stored reaches NB-1 or NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
+ CALL DGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
+ $ VR, LDVR,
+ $ WORK( 1 + (IV)*N ), N,
+ $ ZERO,
+ $ WORK( 1 + (NB+IV)*N ), N )
+* normalize vectors
+ DO K = IV, NB
+ IF( ISCOMPLEX(K).EQ.0 ) THEN
+* real eigenvector
+ II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
+ ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
+* first eigenvector of conjugate pair
+ EMAX = ZERO
+ DO II = 1, N
+ EMAX = MAX( EMAX,
+ $ ABS( WORK( II + (NB+K )*N ) )+
+ $ ABS( WORK( II + (NB+K+1)*N ) ) )
+ END DO
+ REMAX = ONE / EMAX
+* else if ISCOMPLEX(K).EQ.-1
+* second eigenvector of conjugate pair
+* reuse same REMAX as previous K
+ END IF
+ CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL DLACPY( 'F', N, NB-IV+1,
+ $ WORK( 1 + (NB+IV)*N ), N,
+ $ VR( 1, KI2 ), LDVR )
+ IV = NB
+ ELSE
+ IV = IV - 1
+ END IF
+ END IF ! blocked back-transform
+*
+ IS = IS - 1
+ IF( IP.NE.0 )
+ $ IS = IS - 1
+ 140 CONTINUE
+ END IF
+
+ IF( LEFTV ) THEN
+*
+* ============================================================
+* Compute left eigenvectors.
+*
+* IV is index of column in current block.
+* For complex left vector, uses IV for real part and IV+1 for complex part.
+* Non-blocked version always uses IV=1;
+* blocked version starts with IV=1, goes up to NB-1 or NB.
+* (Note the "0-th" column is used for 1-norms computed above.)
+ IV = 1
+ IP = 0
+ IS = 1
+ DO 260 KI = 1, N
+ IF( IP.EQ.1 ) THEN
+* previous iteration (ki-1) was first of conjugate pair,
+* so this ki is second of conjugate pair; skip to end of loop
+ IP = -1
+ GO TO 260
+ ELSE IF( KI.EQ.N ) THEN
+* last column, so this ki must be real eigenvalue
+ IP = 0
+ ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
+* zero on sub-diagonal, so this ki is real eigenvalue
+ IP = 0
+ ELSE
+* non-zero on sub-diagonal, so this ki is first of conjugate pair
+ IP = 1
+ END IF
+*
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 260
+ END IF
+*
+* Compute the KI-th eigenvalue (WR,WI).
+*
+ WR = T( KI, KI )
+ WI = ZERO
+ IF( IP.NE.0 )
+ $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
+ $ SQRT( ABS( T( KI+1, KI ) ) )
+ SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
+*
+ IF( IP.EQ.0 ) THEN
+*
+* --------------------------------------------------------
+* Real left eigenvector
+*
+ WORK( KI + IV*N ) = ONE
+*
+* Form right-hand side.
+*
+ DO 160 K = KI + 1, N
+ WORK( K + IV*N ) = -T( KI, K )
+ 160 CONTINUE
+*
+* Solve transposed quasi-triangular system:
+* [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
+*
+ VMAX = ONE
+ VCRIT = BIGNUM
+*
+ JNXT = KI + 1
+ DO 170 J = KI + 1, N
+ IF( J.LT.JNXT )
+ $ GO TO 170
+ J1 = J
+ J2 = J
+ JNXT = J + 1
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).NE.ZERO ) THEN
+ J2 = J + 1
+ JNXT = J + 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side.
+*
+ IF( WORK( J ).GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+IV*N ) = WORK( J+IV*N ) -
+ $ DDOT( J-KI-1, T( KI+1, J ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+* Solve [ T(J,J) - WR ]**T * X = WORK
+*
+ CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
+ WORK( J+IV*N ) = X( 1, 1 )
+ VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side.
+*
+ BETA = MAX( WORK( J ), WORK( J+1 ) )
+ IF( BETA.GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL DSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+IV*N ) = WORK( J+IV*N ) -
+ $ DDOT( J-KI-1, T( KI+1, J ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+ WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
+ $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+* Solve
+* [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
+* [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
+*
+ CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
+ WORK( J +IV*N ) = X( 1, 1 )
+ WORK( J+1+IV*N ) = X( 2, 1 )
+*
+ VMAX = MAX( ABS( WORK( J +IV*N ) ),
+ $ ABS( WORK( J+1+IV*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ END IF
+ 170 CONTINUE
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL DCOPY( N-KI+1, WORK( KI + IV*N ), 1,
+ $ VL( KI, IS ), 1 )
+*
+ II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
+ REMAX = ONE / ABS( VL( II, IS ) )
+ CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+*
+ DO 180 K = 1, KI - 1
+ VL( K, IS ) = ZERO
+ 180 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N )
+ $ CALL DGEMV( 'N', N, N-KI, ONE,
+ $ VL( 1, KI+1 ), LDVL,
+ $ WORK( KI+1 + IV*N ), 1,
+ $ WORK( KI + IV*N ), VL( 1, KI ), 1 )
+*
+ II = IDAMAX( N, VL( 1, KI ), 1 )
+ REMAX = ONE / ABS( VL( II, KI ) )
+ CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + IV*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+* back-transform and normalization is done below
+ END IF
+ ELSE
+*
+* --------------------------------------------------------
+* Complex left eigenvector.
+*
+* Initial solve:
+* [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
+* [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
+*
+ IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
+ WORK( KI + (IV )*N ) = WI / T( KI, KI+1 )
+ WORK( KI+1 + (IV+1)*N ) = ONE
+ ELSE
+ WORK( KI + (IV )*N ) = ONE
+ WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
+ END IF
+ WORK( KI+1 + (IV )*N ) = ZERO
+ WORK( KI + (IV+1)*N ) = ZERO
+*
+* Form right-hand side.
+*
+ DO 190 K = KI + 2, N
+ WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(KI, K)
+ WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
+ 190 CONTINUE
+*
+* Solve transposed quasi-triangular system:
+* [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
+*
+ VMAX = ONE
+ VCRIT = BIGNUM
+*
+ JNXT = KI + 2
+ DO 200 J = KI + 2, N
+ IF( J.LT.JNXT )
+ $ GO TO 200
+ J1 = J
+ J2 = J
+ JNXT = J + 1
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).NE.ZERO ) THEN
+ J2 = J + 1
+ JNXT = J + 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+* Scale if necessary to avoid overflow when
+* forming the right-hand side elements.
+*
+ IF( WORK( J ).GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
+ CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+(IV )*N ) = WORK( J+(IV)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+ WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+* Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
+*
+ CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ -WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
+ CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
+ END IF
+ WORK( J+(IV )*N ) = X( 1, 1 )
+ WORK( J+(IV+1)*N ) = X( 1, 2 )
+ VMAX = MAX( ABS( WORK( J+(IV )*N ) ),
+ $ ABS( WORK( J+(IV+1)*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side elements.
+*
+ BETA = MAX( WORK( J ), WORK( J+1 ) )
+ IF( BETA.GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL DSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
+ CALL DSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J +(IV )*N ) = WORK( J+(IV)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+*
+ WORK( J +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+ WORK( J+1+(IV )*N ) = WORK( J+1+(IV)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+*
+ WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
+ $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+* Solve 2-by-2 complex linear equation
+* [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
+* [ (T(j+1,j) T(j+1,j+1)) ]
+*
+ CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ -WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
+ CALL DSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
+ END IF
+ WORK( J +(IV )*N ) = X( 1, 1 )
+ WORK( J +(IV+1)*N ) = X( 1, 2 )
+ WORK( J+1+(IV )*N ) = X( 2, 1 )
+ WORK( J+1+(IV+1)*N ) = X( 2, 2 )
+ VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
+ $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
+ $ VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ END IF
+ 200 CONTINUE
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL DCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
+ $ VL( KI, IS ), 1 )
+ CALL DCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
+ $ VL( KI, IS+1 ), 1 )
+*
+ EMAX = ZERO
+ DO 220 K = KI, N
+ EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
+ $ ABS( VL( K, IS+1 ) ) )
+ 220 CONTINUE
+ REMAX = ONE / EMAX
+ CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+ CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
+*
+ DO 230 K = 1, KI - 1
+ VL( K, IS ) = ZERO
+ VL( K, IS+1 ) = ZERO
+ 230 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N-1 ) THEN
+ CALL DGEMV( 'N', N, N-KI-1, ONE,
+ $ VL( 1, KI+2 ), LDVL,
+ $ WORK( KI+2 + (IV)*N ), 1,
+ $ WORK( KI + (IV)*N ),
+ $ VL( 1, KI ), 1 )
+ CALL DGEMV( 'N', N, N-KI-1, ONE,
+ $ VL( 1, KI+2 ), LDVL,
+ $ WORK( KI+2 + (IV+1)*N ), 1,
+ $ WORK( KI+1 + (IV+1)*N ),
+ $ VL( 1, KI+1 ), 1 )
+ ELSE
+ CALL DSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
+ CALL DSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
+ END IF
+*
+ EMAX = ZERO
+ DO 240 K = 1, N
+ EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
+ $ ABS( VL( K, KI+1 ) ) )
+ 240 CONTINUE
+ REMAX = ONE / EMAX
+ CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
+ CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + (IV )*N ) = ZERO
+ WORK( K + (IV+1)*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+ ISCOMPLEX( IV+1 ) = -IP
+ IV = IV + 1
+* back-transform and normalization is done below
+ END IF
+ END IF
+
+ IF( NB.GT.1 ) THEN
+* --------------------------------------------------------
+* Blocked version of back-transform
+* For complex case, KI2 includes both vectors (KI and KI+1)
+ IF( IP.EQ.0 ) THEN
+ KI2 = KI
+ ELSE
+ KI2 = KI + 1
+ END IF
+
+* Columns 1:IV of work are valid vectors.
+* When the number of vectors stored reaches NB-1 or NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
+ CALL DGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
+ $ VL( 1, KI2-IV+1 ), LDVL,
+ $ WORK( KI2-IV+1 + (1)*N ), N,
+ $ ZERO,
+ $ WORK( 1 + (NB+1)*N ), N )
+* normalize vectors
+ DO K = 1, IV
+ IF( ISCOMPLEX(K).EQ.0) THEN
+* real eigenvector
+ II = IDAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
+ ELSE IF( ISCOMPLEX(K).EQ.1) THEN
+* first eigenvector of conjugate pair
+ EMAX = ZERO
+ DO II = 1, N
+ EMAX = MAX( EMAX,
+ $ ABS( WORK( II + (NB+K )*N ) )+
+ $ ABS( WORK( II + (NB+K+1)*N ) ) )
+ END DO
+ REMAX = ONE / EMAX
+* else if ISCOMPLEX(K).EQ.-1
+* second eigenvector of conjugate pair
+* reuse same REMAX as previous K
+ END IF
+ CALL DSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL DLACPY( 'F', N, IV,
+ $ WORK( 1 + (NB+1)*N ), N,
+ $ VL( 1, KI2-IV+1 ), LDVL )
+ IV = 1
+ ELSE
+ IV = IV + 1
+ END IF
+ END IF ! blocked back-transform
+*
+ IS = IS + 1
+ IF( IP.NE.0 )
+ $ IS = IS + 1
+ 260 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of DTREVC3
+*
+ END
diff --git a/SRC/ilaenv.f b/SRC/ilaenv.f
index 89a4468f..a8daf00a 100644
--- a/SRC/ilaenv.f
+++ b/SRC/ilaenv.f
@@ -397,6 +397,12 @@
ELSE
NB = 64
END IF
+ ELSE IF ( C3.EQ.'EVC' ) THEN
+ IF( SNAME ) THEN
+ NB = 64
+ ELSE
+ NB = 64
+ END IF
END IF
ELSE IF( C2.EQ.'LA' ) THEN
IF( C3.EQ.'UUM' ) THEN
diff --git a/SRC/sgeev.f b/SRC/sgeev.f
index 667de0af..1187f5c3 100644
--- a/SRC/sgeev.f
+++ b/SRC/sgeev.f
@@ -188,6 +188,7 @@
* =====================================================================
SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
$ LDVR, WORK, LWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -213,7 +214,7 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
CHARACTER SIDE
INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
- $ MAXWRK, MINWRK, NOUT
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
$ SN
* ..
@@ -223,7 +224,7 @@
* ..
* .. External Subroutines ..
EXTERNAL SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY,
- $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC,
+ $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC3,
$ XERBLA
* ..
* .. External Functions ..
@@ -279,24 +280,34 @@
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
$ 'SORGHR', ' ', N, 1, N, -1 ) )
CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
+ CALL STREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR, N, NOUT,
+ $ WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE IF( WANTVR ) THEN
MINWRK = 4*N
MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
$ 'SORGHR', ' ', N, 1, N, -1 ) )
CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
+ CALL STREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR, N, NOUT,
+ $ WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
MAXWRK = MAX( MAXWRK, 4*N )
ELSE
MINWRK = 3*N
CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR,
- $ WORK, -1, INFO )
- HSWORK = WORK( 1 )
+ $ WORK, -1, INFO )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK )
END IF
MAXWRK = MAX( MAXWRK, MINWRK )
@@ -426,10 +437,10 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (Workspace: need 4*N)
+* (Workspace: need 4*N, prefer N + N + 2*N*NB)
*
- CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), IERR )
+ CALL STREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
END IF
*
IF( WANTVL ) THEN
diff --git a/SRC/sgeevx.f b/SRC/sgeevx.f
index 93f63b0b..eff3a9f4 100644
--- a/SRC/sgeevx.f
+++ b/SRC/sgeevx.f
@@ -302,6 +302,7 @@
SUBROUTINE SGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
$ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
$ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.2) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -330,8 +331,8 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
$ WNTSNN, WNTSNV
CHARACTER JOB, SIDE
- INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
- $ MINWRK, NOUT
+ INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
$ SN
* ..
@@ -341,7 +342,7 @@
* ..
* .. External Subroutines ..
EXTERNAL SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY,
- $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC,
+ $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC3,
$ STRSNA, XERBLA
* ..
* .. External Functions ..
@@ -366,8 +367,9 @@
WNTSNE = LSAME( SENSE, 'E' )
WNTSNV = LSAME( SENSE, 'V' )
WNTSNB = LSAME( SENSE, 'B' )
- IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR.
- $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN
+ IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' )
+ $ .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) )
+ $ THEN
INFO = -1
ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
INFO = -2
@@ -405,9 +407,19 @@
MAXWRK = N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 )
*
IF( WANTVL ) THEN
+ CALL STREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL,
$ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
+ CALL STREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR,
$ WORK, -1, INFO )
ELSE
@@ -419,7 +431,7 @@
$ LDVR, WORK, -1, INFO )
END IF
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
*
IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
MINWRK = 2*N
@@ -579,10 +591,10 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (Workspace: need 3*N)
+* (Workspace: need 3*N, prefer N + 2*N*NB)
*
- CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), IERR )
+ CALL STREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR )
END IF
*
* Compute condition numbers if desired
diff --git a/SRC/strevc3.f b/SRC/strevc3.f
new file mode 100644
index 00000000..95ac0f6d
--- /dev/null
+++ b/SRC/strevc3.f
@@ -0,0 +1,1303 @@
+*> \brief \b STREVC3
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download STREVC3 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strevc3.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strevc3.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strevc3.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE STREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+* VR, LDVR, MM, M, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER HOWMNY, SIDE
+* INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+* $ WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> STREVC3 computes some or all of the right and/or left eigenvectors of
+*> a real upper quasi-triangular matrix T.
+*> Matrices of this type are produced by the Schur factorization of
+*> a real general matrix: A = Q*T*Q**T, as computed by SHSEQR.
+*>
+*> The right eigenvector x and the left eigenvector y of T corresponding
+*> to an eigenvalue w are defined by:
+*>
+*> T*x = w*x, (y**T)*T = w*(y**T)
+*>
+*> where y**T denotes the transpose of the vector y.
+*> The eigenvalues are not input to this routine, but are read directly
+*> from the diagonal blocks 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 orthogonal 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*> SIDE is CHARACTER*1
+*> = 'R': compute right eigenvectors only;
+*> = 'L': compute left eigenvectors only;
+*> = 'B': compute both right and left eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] HOWMNY
+*> \verbatim
+*> HOWMNY is CHARACTER*1
+*> = 'A': compute all right and/or left eigenvectors;
+*> = 'B': compute all right and/or left eigenvectors,
+*> backtransformed by the matrices in VR and/or VL;
+*> = 'S': compute selected right and/or left eigenvectors,
+*> as indicated by the logical array SELECT.
+*> \endverbatim
+*>
+*> \param[in,out] SELECT
+*> \verbatim
+*> SELECT is LOGICAL array, dimension (N)
+*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
+*> computed.
+*> If w(j) is a real eigenvalue, the corresponding real
+*> eigenvector is computed if SELECT(j) is .TRUE..
+*> If w(j) and w(j+1) are the real and imaginary parts of a
+*> complex eigenvalue, the corresponding complex eigenvector is
+*> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
+*> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
+*> .FALSE..
+*> Not referenced if HOWMNY = 'A' or 'B'.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix T. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*> T is REAL array, dimension (LDT,N)
+*> The upper quasi-triangular matrix T in Schur canonical form.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*> LDT is INTEGER
+*> The leading dimension of the array T. LDT >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] VL
+*> \verbatim
+*> VL is REAL 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 orthogonal matrix Q
+*> of Schur vectors returned by SHSEQR).
+*> 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.
+*> A complex eigenvector corresponding to a complex eigenvalue
+*> is stored in two consecutive columns, the first holding the
+*> real part, and the second the imaginary part.
+*> Not referenced if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the array VL.
+*> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[in,out] VR
+*> \verbatim
+*> VR is REAL 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 orthogonal matrix Q
+*> of Schur vectors returned by SHSEQR).
+*> 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.
+*> A complex eigenvector corresponding to a complex eigenvalue
+*> is stored in two consecutive columns, the first holding the
+*> real part and the second the imaginary part.
+*> Not referenced if SIDE = 'L'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the array VR.
+*> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[in] MM
+*> \verbatim
+*> MM is INTEGER
+*> The number of columns in the arrays VL and/or VR. MM >= M.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> 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 real eigenvector occupies one column and each
+*> selected complex eigenvector occupies two columns.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is REAL array, dimension (MAX(1,LWORK))
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of array WORK. LWORK >= max(1,3*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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+* @generated from dtrevc3.f, fortran d -> s, Tue Apr 19 01:47:44 2016
+*
+*> \ingroup realOTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> 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|.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE STREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+ $ VR, LDVR, MM, M, WORK, LWORK, INFO )
+ IMPLICIT NONE
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER HOWMNY, SIDE
+ INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+ LOGICAL SELECT( * )
+ REAL T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+ $ WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ REAL ZERO, ONE
+ PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
+ INTEGER NBMIN, NBMAX
+ PARAMETER ( NBMIN = 8, NBMAX = 128 )
+* ..
+* .. Local Scalars ..
+ LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
+ $ RIGHTV, SOMEV
+ INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
+ $ IV, MAXWRK, NB, KI2
+ REAL BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
+ $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
+ $ XNORM
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ISAMAX, ILAENV
+ REAL SDOT, SLAMCH
+ EXTERNAL LSAME, ISAMAX, ILAENV, SDOT, SLAMCH
+* ..
+* .. External Subroutines ..
+ EXTERNAL SAXPY, SCOPY, SGEMV, SLALN2, SSCAL, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, SQRT
+* ..
+* .. Local Arrays ..
+ REAL X( 2, 2 )
+ INTEGER ISCOMPLEX( NBMAX )
+* ..
+* .. Executable Statements ..
+*
+* Decode and test the input parameters
+*
+ BOTHV = LSAME( SIDE, 'B' )
+ RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
+ LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
+*
+ ALLV = LSAME( HOWMNY, 'A' )
+ OVER = LSAME( HOWMNY, 'B' )
+ SOMEV = LSAME( HOWMNY, 'S' )
+*
+ INFO = 0
+ NB = ILAENV( 1, 'STREVC', SIDE // HOWMNY, N, -1, -1, -1 )
+ MAXWRK = N + 2*N*NB
+ WORK(1) = MAXWRK
+ LQUERY = ( LWORK.EQ.-1 )
+ IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
+ INFO = -1
+ ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
+ INFO = -8
+ ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
+ INFO = -10
+ ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -14
+ ELSE
+*
+* Set M to the number of columns required to store the selected
+* eigenvectors, standardize the array SELECT if necessary, and
+* test MM.
+*
+ IF( SOMEV ) THEN
+ M = 0
+ PAIR = .FALSE.
+ DO 10 J = 1, N
+ IF( PAIR ) THEN
+ PAIR = .FALSE.
+ SELECT( J ) = .FALSE.
+ ELSE
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).EQ.ZERO ) THEN
+ IF( SELECT( J ) )
+ $ M = M + 1
+ ELSE
+ PAIR = .TRUE.
+ IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
+ SELECT( J ) = .TRUE.
+ M = M + 2
+ END IF
+ END IF
+ ELSE
+ IF( SELECT( N ) )
+ $ M = M + 1
+ END IF
+ END IF
+ 10 CONTINUE
+ ELSE
+ M = N
+ END IF
+*
+ IF( MM.LT.M ) THEN
+ INFO = -11
+ END IF
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'STREVC3', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Use blocked version of back-transformation if sufficient workspace.
+* Zero-out the workspace to avoid potential NaN propagation.
+*
+ IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
+ NB = (LWORK - N) / (2*N)
+ NB = MIN( NB, NBMAX )
+ CALL SLASET( 'F', N, 1+2*NB, ZERO, ZERO, WORK, N )
+ ELSE
+ NB = 1
+ END IF
+*
+* Set the constants to control overflow.
+*
+ UNFL = SLAMCH( 'Safe minimum' )
+ OVFL = ONE / UNFL
+ CALL SLABAD( UNFL, OVFL )
+ ULP = SLAMCH( 'Precision' )
+ SMLNUM = UNFL*( N / ULP )
+ BIGNUM = ( ONE-ULP ) / SMLNUM
+*
+* Compute 1-norm of each column of strictly upper triangular
+* part of T to control overflow in triangular solver.
+*
+ WORK( 1 ) = ZERO
+ DO 30 J = 2, N
+ WORK( J ) = ZERO
+ DO 20 I = 1, J - 1
+ WORK( J ) = WORK( J ) + ABS( T( I, J ) )
+ 20 CONTINUE
+ 30 CONTINUE
+*
+* Index IP is used to specify the real or complex eigenvalue:
+* IP = 0, real eigenvalue,
+* 1, first of conjugate complex pair: (wr,wi)
+* -1, second of conjugate complex pair: (wr,wi)
+* ISCOMPLEX array stores IP for each column in current block.
+*
+ IF( RIGHTV ) THEN
+*
+* ============================================================
+* Compute right eigenvectors.
+*
+* IV is index of column in current block.
+* For complex right vector, uses IV-1 for real part and IV for complex part.
+* Non-blocked version always uses IV=2;
+* blocked version starts with IV=NB, goes down to 1 or 2.
+* (Note the "0-th" column is used for 1-norms computed above.)
+ IV = 2
+ IF( NB.GT.2 ) THEN
+ IV = NB
+ END IF
+
+ IP = 0
+ IS = M
+ DO 140 KI = N, 1, -1
+ IF( IP.EQ.-1 ) THEN
+* previous iteration (ki+1) was second of conjugate pair,
+* so this ki is first of conjugate pair; skip to end of loop
+ IP = 1
+ GO TO 140
+ ELSE IF( KI.EQ.1 ) THEN
+* last column, so this ki must be real eigenvalue
+ IP = 0
+ ELSE IF( T( KI, KI-1 ).EQ.ZERO ) THEN
+* zero on sub-diagonal, so this ki is real eigenvalue
+ IP = 0
+ ELSE
+* non-zero on sub-diagonal, so this ki is second of conjugate pair
+ IP = -1
+ END IF
+
+ IF( SOMEV ) THEN
+ IF( IP.EQ.0 ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 140
+ ELSE
+ IF( .NOT.SELECT( KI-1 ) )
+ $ GO TO 140
+ END IF
+ END IF
+*
+* Compute the KI-th eigenvalue (WR,WI).
+*
+ WR = T( KI, KI )
+ WI = ZERO
+ IF( IP.NE.0 )
+ $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
+ $ SQRT( ABS( T( KI-1, KI ) ) )
+ SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
+*
+ IF( IP.EQ.0 ) THEN
+*
+* --------------------------------------------------------
+* Real right eigenvector
+*
+ WORK( KI + IV*N ) = ONE
+*
+* Form right-hand side.
+*
+ DO 50 K = 1, KI - 1
+ WORK( K + IV*N ) = -T( K, KI )
+ 50 CONTINUE
+*
+* Solve upper quasi-triangular system:
+* [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
+*
+ JNXT = KI - 1
+ DO 60 J = KI - 1, 1, -1
+ IF( J.GT.JNXT )
+ $ GO TO 60
+ J1 = J
+ J2 = J
+ JNXT = J - 1
+ IF( J.GT.1 ) THEN
+ IF( T( J, J-1 ).NE.ZERO ) THEN
+ J1 = J - 1
+ JNXT = J - 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+ CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale X(1,1) to avoid overflow when updating
+* the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL SSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
+ WORK( J+IV*N ) = X( 1, 1 )
+*
+* Update right-hand side
+*
+ CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
+ $ WORK( 1+IV*N ), 1 )
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+ CALL SLALN2( .FALSE., 2, 1, SMIN, ONE,
+ $ T( J-1, J-1 ), LDT, ONE, ONE,
+ $ WORK( J-1+IV*N ), N, WR, ZERO, X, 2,
+ $ SCALE, XNORM, IERR )
+*
+* Scale X(1,1) and X(2,1) to avoid overflow when
+* updating the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ BETA = MAX( WORK( J-1 ), WORK( J ) )
+ IF( BETA.GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ X( 2, 1 ) = X( 2, 1 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL SSCAL( KI, SCALE, WORK( 1+IV*N ), 1 )
+ WORK( J-1+IV*N ) = X( 1, 1 )
+ WORK( J +IV*N ) = X( 2, 1 )
+*
+* Update right-hand side
+*
+ CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
+ $ WORK( 1+IV*N ), 1 )
+ CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
+ $ WORK( 1+IV*N ), 1 )
+ END IF
+ 60 CONTINUE
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL SCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
+*
+ II = ISAMAX( KI, VR( 1, IS ), 1 )
+ REMAX = ONE / ABS( VR( II, IS ) )
+ CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 70 K = KI + 1, N
+ VR( K, IS ) = ZERO
+ 70 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.1 )
+ $ CALL SGEMV( 'N', N, KI-1, ONE, VR, LDVR,
+ $ WORK( 1 + IV*N ), 1, WORK( KI + IV*N ),
+ $ VR( 1, KI ), 1 )
+*
+ II = ISAMAX( N, VR( 1, KI ), 1 )
+ REMAX = ONE / ABS( VR( II, KI ) )
+ CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + IV*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+* back-transform and normalization is done below
+ END IF
+ ELSE
+*
+* --------------------------------------------------------
+* Complex right eigenvector.
+*
+* Initial solve
+* [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
+* [ ( T(KI, KI-1) T(KI, KI) ) ]
+*
+ IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
+ WORK( KI-1 + (IV-1)*N ) = ONE
+ WORK( KI + (IV )*N ) = WI / T( KI-1, KI )
+ ELSE
+ WORK( KI-1 + (IV-1)*N ) = -WI / T( KI, KI-1 )
+ WORK( KI + (IV )*N ) = ONE
+ END IF
+ WORK( KI + (IV-1)*N ) = ZERO
+ WORK( KI-1 + (IV )*N ) = ZERO
+*
+* Form right-hand side.
+*
+ DO 80 K = 1, KI - 2
+ WORK( K+(IV-1)*N ) = -WORK( KI-1+(IV-1)*N )*T(K,KI-1)
+ WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(K,KI )
+ 80 CONTINUE
+*
+* Solve upper quasi-triangular system:
+* [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
+*
+ JNXT = KI - 2
+ DO 90 J = KI - 2, 1, -1
+ IF( J.GT.JNXT )
+ $ GO TO 90
+ J1 = J
+ J2 = J
+ JNXT = J - 1
+ IF( J.GT.1 ) THEN
+ IF( T( J, J-1 ).NE.ZERO ) THEN
+ J1 = J - 1
+ JNXT = J - 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+ CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+(IV-1)*N ), N,
+ $ WR, WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale X(1,1) and X(1,2) to avoid overflow when
+* updating the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
+ X( 1, 1 ) = X( 1, 1 ) / XNORM
+ X( 1, 2 ) = X( 1, 2 ) / XNORM
+ SCALE = SCALE / XNORM
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL SSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
+ CALL SSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
+ END IF
+ WORK( J+(IV-1)*N ) = X( 1, 1 )
+ WORK( J+(IV )*N ) = X( 1, 2 )
+*
+* Update the right-hand side
+*
+ CALL SAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL SAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+ CALL SLALN2( .FALSE., 2, 2, SMIN, ONE,
+ $ T( J-1, J-1 ), LDT, ONE, ONE,
+ $ WORK( J-1+(IV-1)*N ), N, WR, WI, X, 2,
+ $ SCALE, XNORM, IERR )
+*
+* Scale X to avoid overflow when updating
+* the right-hand side.
+*
+ IF( XNORM.GT.ONE ) THEN
+ BETA = MAX( WORK( J-1 ), WORK( J ) )
+ IF( BETA.GT.BIGNUM / XNORM ) THEN
+ REC = ONE / XNORM
+ X( 1, 1 ) = X( 1, 1 )*REC
+ X( 1, 2 ) = X( 1, 2 )*REC
+ X( 2, 1 ) = X( 2, 1 )*REC
+ X( 2, 2 ) = X( 2, 2 )*REC
+ SCALE = SCALE*REC
+ END IF
+ END IF
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL SSCAL( KI, SCALE, WORK( 1+(IV-1)*N ), 1 )
+ CALL SSCAL( KI, SCALE, WORK( 1+(IV )*N ), 1 )
+ END IF
+ WORK( J-1+(IV-1)*N ) = X( 1, 1 )
+ WORK( J +(IV-1)*N ) = X( 2, 1 )
+ WORK( J-1+(IV )*N ) = X( 1, 2 )
+ WORK( J +(IV )*N ) = X( 2, 2 )
+*
+* Update the right-hand side
+*
+ CALL SAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL SAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
+ $ WORK( 1+(IV-1)*N ), 1 )
+ CALL SAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+ CALL SAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
+ $ WORK( 1+(IV )*N ), 1 )
+ END IF
+ 90 CONTINUE
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL SCOPY( KI, WORK( 1+(IV-1)*N ), 1, VR(1,IS-1), 1 )
+ CALL SCOPY( KI, WORK( 1+(IV )*N ), 1, VR(1,IS ), 1 )
+*
+ EMAX = ZERO
+ DO 100 K = 1, KI
+ EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
+ $ ABS( VR( K, IS ) ) )
+ 100 CONTINUE
+ REMAX = ONE / EMAX
+ CALL SSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
+ CALL SSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 110 K = KI + 1, N
+ VR( K, IS-1 ) = ZERO
+ VR( K, IS ) = ZERO
+ 110 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.2 ) THEN
+ CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
+ $ WORK( 1 + (IV-1)*N ), 1,
+ $ WORK( KI-1 + (IV-1)*N ), VR(1,KI-1), 1)
+ CALL SGEMV( 'N', N, KI-2, ONE, VR, LDVR,
+ $ WORK( 1 + (IV)*N ), 1,
+ $ WORK( KI + (IV)*N ), VR( 1, KI ), 1 )
+ ELSE
+ CALL SSCAL( N, WORK(KI-1+(IV-1)*N), VR(1,KI-1), 1)
+ CALL SSCAL( N, WORK(KI +(IV )*N), VR(1,KI ), 1)
+ END IF
+*
+ EMAX = ZERO
+ DO 120 K = 1, N
+ EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
+ $ ABS( VR( K, KI ) ) )
+ 120 CONTINUE
+ REMAX = ONE / EMAX
+ CALL SSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
+ CALL SSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + (IV-1)*N ) = ZERO
+ WORK( K + (IV )*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV-1 ) = -IP
+ ISCOMPLEX( IV ) = IP
+ IV = IV - 1
+* back-transform and normalization is done below
+ END IF
+ END IF
+
+ IF( NB.GT.1 ) THEN
+* --------------------------------------------------------
+* Blocked version of back-transform
+* For complex case, KI2 includes both vectors (KI-1 and KI)
+ IF( IP.EQ.0 ) THEN
+ KI2 = KI
+ ELSE
+ KI2 = KI - 1
+ END IF
+
+* Columns IV:NB of work are valid vectors.
+* When the number of vectors stored reaches NB-1 or NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.LE.2) .OR. (KI2.EQ.1) ) THEN
+ CALL SGEMM( 'N', 'N', N, NB-IV+1, KI2+NB-IV, ONE,
+ $ VR, LDVR,
+ $ WORK( 1 + (IV)*N ), N,
+ $ ZERO,
+ $ WORK( 1 + (NB+IV)*N ), N )
+* normalize vectors
+ DO K = IV, NB
+ IF( ISCOMPLEX(K).EQ.0 ) THEN
+* real eigenvector
+ II = ISAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
+ ELSE IF( ISCOMPLEX(K).EQ.1 ) THEN
+* first eigenvector of conjugate pair
+ EMAX = ZERO
+ DO II = 1, N
+ EMAX = MAX( EMAX,
+ $ ABS( WORK( II + (NB+K )*N ) )+
+ $ ABS( WORK( II + (NB+K+1)*N ) ) )
+ END DO
+ REMAX = ONE / EMAX
+* else if ISCOMPLEX(K).EQ.-1
+* second eigenvector of conjugate pair
+* reuse same REMAX as previous K
+ END IF
+ CALL SSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL SLACPY( 'F', N, NB-IV+1,
+ $ WORK( 1 + (NB+IV)*N ), N,
+ $ VR( 1, KI2 ), LDVR )
+ IV = NB
+ ELSE
+ IV = IV - 1
+ END IF
+ END IF ! blocked back-transform
+*
+ IS = IS - 1
+ IF( IP.NE.0 )
+ $ IS = IS - 1
+ 140 CONTINUE
+ END IF
+
+ IF( LEFTV ) THEN
+*
+* ============================================================
+* Compute left eigenvectors.
+*
+* IV is index of column in current block.
+* For complex left vector, uses IV for real part and IV+1 for complex part.
+* Non-blocked version always uses IV=1;
+* blocked version starts with IV=1, goes up to NB-1 or NB.
+* (Note the "0-th" column is used for 1-norms computed above.)
+ IV = 1
+ IP = 0
+ IS = 1
+ DO 260 KI = 1, N
+ IF( IP.EQ.1 ) THEN
+* previous iteration (ki-1) was first of conjugate pair,
+* so this ki is second of conjugate pair; skip to end of loop
+ IP = -1
+ GO TO 260
+ ELSE IF( KI.EQ.N ) THEN
+* last column, so this ki must be real eigenvalue
+ IP = 0
+ ELSE IF( T( KI+1, KI ).EQ.ZERO ) THEN
+* zero on sub-diagonal, so this ki is real eigenvalue
+ IP = 0
+ ELSE
+* non-zero on sub-diagonal, so this ki is first of conjugate pair
+ IP = 1
+ END IF
+*
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 260
+ END IF
+*
+* Compute the KI-th eigenvalue (WR,WI).
+*
+ WR = T( KI, KI )
+ WI = ZERO
+ IF( IP.NE.0 )
+ $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
+ $ SQRT( ABS( T( KI+1, KI ) ) )
+ SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
+*
+ IF( IP.EQ.0 ) THEN
+*
+* --------------------------------------------------------
+* Real left eigenvector
+*
+ WORK( KI + IV*N ) = ONE
+*
+* Form right-hand side.
+*
+ DO 160 K = KI + 1, N
+ WORK( K + IV*N ) = -T( KI, K )
+ 160 CONTINUE
+*
+* Solve transposed quasi-triangular system:
+* [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
+*
+ VMAX = ONE
+ VCRIT = BIGNUM
+*
+ JNXT = KI + 1
+ DO 170 J = KI + 1, N
+ IF( J.LT.JNXT )
+ $ GO TO 170
+ J1 = J
+ J2 = J
+ JNXT = J + 1
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).NE.ZERO ) THEN
+ J2 = J + 1
+ JNXT = J + 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side.
+*
+ IF( WORK( J ).GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL SSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+IV*N ) = WORK( J+IV*N ) -
+ $ SDOT( J-KI-1, T( KI+1, J ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+* Solve [ T(J,J) - WR ]**T * X = WORK
+*
+ CALL SLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
+ WORK( J+IV*N ) = X( 1, 1 )
+ VMAX = MAX( ABS( WORK( J+IV*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side.
+*
+ BETA = MAX( WORK( J ), WORK( J+1 ) )
+ IF( BETA.GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL SSCAL( N-KI+1, REC, WORK( KI+IV*N ), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+IV*N ) = WORK( J+IV*N ) -
+ $ SDOT( J-KI-1, T( KI+1, J ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+ WORK( J+1+IV*N ) = WORK( J+1+IV*N ) -
+ $ SDOT( J-KI-1, T( KI+1, J+1 ), 1,
+ $ WORK( KI+1+IV*N ), 1 )
+*
+* Solve
+* [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
+* [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
+*
+ CALL SLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ ZERO, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE )
+ $ CALL SSCAL( N-KI+1, SCALE, WORK( KI+IV*N ), 1 )
+ WORK( J +IV*N ) = X( 1, 1 )
+ WORK( J+1+IV*N ) = X( 2, 1 )
+*
+ VMAX = MAX( ABS( WORK( J +IV*N ) ),
+ $ ABS( WORK( J+1+IV*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ END IF
+ 170 CONTINUE
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL SCOPY( N-KI+1, WORK( KI + IV*N ), 1,
+ $ VL( KI, IS ), 1 )
+*
+ II = ISAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
+ REMAX = ONE / ABS( VL( II, IS ) )
+ CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+*
+ DO 180 K = 1, KI - 1
+ VL( K, IS ) = ZERO
+ 180 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N )
+ $ CALL SGEMV( 'N', N, N-KI, ONE,
+ $ VL( 1, KI+1 ), LDVL,
+ $ WORK( KI+1 + IV*N ), 1,
+ $ WORK( KI + IV*N ), VL( 1, KI ), 1 )
+*
+ II = ISAMAX( N, VL( 1, KI ), 1 )
+ REMAX = ONE / ABS( VL( II, KI ) )
+ CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + IV*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+* back-transform and normalization is done below
+ END IF
+ ELSE
+*
+* --------------------------------------------------------
+* Complex left eigenvector.
+*
+* Initial solve:
+* [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
+* [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
+*
+ IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
+ WORK( KI + (IV )*N ) = WI / T( KI, KI+1 )
+ WORK( KI+1 + (IV+1)*N ) = ONE
+ ELSE
+ WORK( KI + (IV )*N ) = ONE
+ WORK( KI+1 + (IV+1)*N ) = -WI / T( KI+1, KI )
+ END IF
+ WORK( KI+1 + (IV )*N ) = ZERO
+ WORK( KI + (IV+1)*N ) = ZERO
+*
+* Form right-hand side.
+*
+ DO 190 K = KI + 2, N
+ WORK( K+(IV )*N ) = -WORK( KI +(IV )*N )*T(KI, K)
+ WORK( K+(IV+1)*N ) = -WORK( KI+1+(IV+1)*N )*T(KI+1,K)
+ 190 CONTINUE
+*
+* Solve transposed quasi-triangular system:
+* [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
+*
+ VMAX = ONE
+ VCRIT = BIGNUM
+*
+ JNXT = KI + 2
+ DO 200 J = KI + 2, N
+ IF( J.LT.JNXT )
+ $ GO TO 200
+ J1 = J
+ J2 = J
+ JNXT = J + 1
+ IF( J.LT.N ) THEN
+ IF( T( J+1, J ).NE.ZERO ) THEN
+ J2 = J + 1
+ JNXT = J + 2
+ END IF
+ END IF
+*
+ IF( J1.EQ.J2 ) THEN
+*
+* 1-by-1 diagonal block
+*
+* Scale if necessary to avoid overflow when
+* forming the right-hand side elements.
+*
+ IF( WORK( J ).GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL SSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
+ CALL SSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J+(IV )*N ) = WORK( J+(IV)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+ WORK( J+(IV+1)*N ) = WORK( J+(IV+1)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+* Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
+*
+ CALL SLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ -WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
+ CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
+ END IF
+ WORK( J+(IV )*N ) = X( 1, 1 )
+ WORK( J+(IV+1)*N ) = X( 1, 2 )
+ VMAX = MAX( ABS( WORK( J+(IV )*N ) ),
+ $ ABS( WORK( J+(IV+1)*N ) ), VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ ELSE
+*
+* 2-by-2 diagonal block
+*
+* Scale if necessary to avoid overflow when forming
+* the right-hand side elements.
+*
+ BETA = MAX( WORK( J ), WORK( J+1 ) )
+ IF( BETA.GT.VCRIT ) THEN
+ REC = ONE / VMAX
+ CALL SSCAL( N-KI+1, REC, WORK(KI+(IV )*N), 1 )
+ CALL SSCAL( N-KI+1, REC, WORK(KI+(IV+1)*N), 1 )
+ VMAX = ONE
+ VCRIT = BIGNUM
+ END IF
+*
+ WORK( J +(IV )*N ) = WORK( J+(IV)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+*
+ WORK( J +(IV+1)*N ) = WORK( J+(IV+1)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+ WORK( J+1+(IV )*N ) = WORK( J+1+(IV)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
+ $ WORK( KI+2+(IV)*N ), 1 )
+*
+ WORK( J+1+(IV+1)*N ) = WORK( J+1+(IV+1)*N ) -
+ $ SDOT( J-KI-2, T( KI+2, J+1 ), 1,
+ $ WORK( KI+2+(IV+1)*N ), 1 )
+*
+* Solve 2-by-2 complex linear equation
+* [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
+* [ (T(j+1,j) T(j+1,j+1)) ]
+*
+ CALL SLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
+ $ LDT, ONE, ONE, WORK( J+IV*N ), N, WR,
+ $ -WI, X, 2, SCALE, XNORM, IERR )
+*
+* Scale if necessary
+*
+ IF( SCALE.NE.ONE ) THEN
+ CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV )*N), 1)
+ CALL SSCAL( N-KI+1, SCALE, WORK(KI+(IV+1)*N), 1)
+ END IF
+ WORK( J +(IV )*N ) = X( 1, 1 )
+ WORK( J +(IV+1)*N ) = X( 1, 2 )
+ WORK( J+1+(IV )*N ) = X( 2, 1 )
+ WORK( J+1+(IV+1)*N ) = X( 2, 2 )
+ VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
+ $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ),
+ $ VMAX )
+ VCRIT = BIGNUM / VMAX
+*
+ END IF
+ 200 CONTINUE
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL SCOPY( N-KI+1, WORK( KI + (IV )*N ), 1,
+ $ VL( KI, IS ), 1 )
+ CALL SCOPY( N-KI+1, WORK( KI + (IV+1)*N ), 1,
+ $ VL( KI, IS+1 ), 1 )
+*
+ EMAX = ZERO
+ DO 220 K = KI, N
+ EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
+ $ ABS( VL( K, IS+1 ) ) )
+ 220 CONTINUE
+ REMAX = ONE / EMAX
+ CALL SSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+ CALL SSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
+*
+ DO 230 K = 1, KI - 1
+ VL( K, IS ) = ZERO
+ VL( K, IS+1 ) = ZERO
+ 230 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N-1 ) THEN
+ CALL SGEMV( 'N', N, N-KI-1, ONE,
+ $ VL( 1, KI+2 ), LDVL,
+ $ WORK( KI+2 + (IV)*N ), 1,
+ $ WORK( KI + (IV)*N ),
+ $ VL( 1, KI ), 1 )
+ CALL SGEMV( 'N', N, N-KI-1, ONE,
+ $ VL( 1, KI+2 ), LDVL,
+ $ WORK( KI+2 + (IV+1)*N ), 1,
+ $ WORK( KI+1 + (IV+1)*N ),
+ $ VL( 1, KI+1 ), 1 )
+ ELSE
+ CALL SSCAL( N, WORK(KI+ (IV )*N), VL(1, KI ), 1)
+ CALL SSCAL( N, WORK(KI+1+(IV+1)*N), VL(1, KI+1), 1)
+ END IF
+*
+ EMAX = ZERO
+ DO 240 K = 1, N
+ EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
+ $ ABS( VL( K, KI+1 ) ) )
+ 240 CONTINUE
+ REMAX = ONE / EMAX
+ CALL SSCAL( N, REMAX, VL( 1, KI ), 1 )
+ CALL SSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + (IV )*N ) = ZERO
+ WORK( K + (IV+1)*N ) = ZERO
+ END DO
+ ISCOMPLEX( IV ) = IP
+ ISCOMPLEX( IV+1 ) = -IP
+ IV = IV + 1
+* back-transform and normalization is done below
+ END IF
+ END IF
+
+ IF( NB.GT.1 ) THEN
+* --------------------------------------------------------
+* Blocked version of back-transform
+* For complex case, KI2 includes both vectors (KI and KI+1)
+ IF( IP.EQ.0 ) THEN
+ KI2 = KI
+ ELSE
+ KI2 = KI + 1
+ END IF
+
+* Columns 1:IV of work are valid vectors.
+* When the number of vectors stored reaches NB-1 or NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.GE.NB-1) .OR. (KI2.EQ.N) ) THEN
+ CALL SGEMM( 'N', 'N', N, IV, N-KI2+IV, ONE,
+ $ VL( 1, KI2-IV+1 ), LDVL,
+ $ WORK( KI2-IV+1 + (1)*N ), N,
+ $ ZERO,
+ $ WORK( 1 + (NB+1)*N ), N )
+* normalize vectors
+ DO K = 1, IV
+ IF( ISCOMPLEX(K).EQ.0) THEN
+* real eigenvector
+ II = ISAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / ABS( WORK( II + (NB+K)*N ) )
+ ELSE IF( ISCOMPLEX(K).EQ.1) THEN
+* first eigenvector of conjugate pair
+ EMAX = ZERO
+ DO II = 1, N
+ EMAX = MAX( EMAX,
+ $ ABS( WORK( II + (NB+K )*N ) )+
+ $ ABS( WORK( II + (NB+K+1)*N ) ) )
+ END DO
+ REMAX = ONE / EMAX
+* else if ISCOMPLEX(K).EQ.-1
+* second eigenvector of conjugate pair
+* reuse same REMAX as previous K
+ END IF
+ CALL SSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL SLACPY( 'F', N, IV,
+ $ WORK( 1 + (NB+1)*N ), N,
+ $ VL( 1, KI2-IV+1 ), LDVL )
+ IV = 1
+ ELSE
+ IV = IV + 1
+ END IF
+ END IF ! blocked back-transform
+*
+ IS = IS + 1
+ IF( IP.NE.0 )
+ $ IS = IS + 1
+ 260 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of STREVC3
+*
+ END
diff --git a/SRC/zgeev.f b/SRC/zgeev.f
index a518b4cd..caed1818 100644
--- a/SRC/zgeev.f
+++ b/SRC/zgeev.f
@@ -176,6 +176,7 @@
* =====================================================================
SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
$ WORK, LWORK, RWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -202,7 +203,7 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
CHARACTER SIDE
INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
- $ IWRK, K, MAXWRK, MINWRK, NOUT
+ $ IWRK, K, LWORK_TREVC, MAXWRK, MINWRK, NOUT
DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
COMPLEX*16 TMP
* ..
@@ -212,7 +213,7 @@
* ..
* .. External Subroutines ..
EXTERNAL DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
- $ ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
+ $ ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC3, ZUNGHR
* ..
* .. External Functions ..
LOGICAL LSAME
@@ -266,18 +267,28 @@
IF( WANTVL ) THEN
MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
$ ' ', N, 1, N, -1 ) )
+ CALL ZTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
$ ' ', N, 1, N, -1 ) )
+ CALL ZTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, N + LWORK_TREVC )
CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
ELSE
CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
- $ WORK, -1, INFO )
+ $ WORK, -1, INFO )
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
END IF
WORK( 1 ) = MAXWRK
@@ -412,12 +423,13 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (CWorkspace: need 2*N)
+* (CWorkspace: need 2*N, prefer N + 2*N*NB)
* (RWorkspace: need 2*N)
*
IRWORK = IBAL + N
- CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
+ CALL ZTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
+ $ RWORK( IRWORK ), N, IERR )
END IF
*
IF( WANTVL ) THEN
diff --git a/SRC/zgeevx.f b/SRC/zgeevx.f
index e68165e2..cb750650 100644
--- a/SRC/zgeevx.f
+++ b/SRC/zgeevx.f
@@ -284,6 +284,7 @@
SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
$ LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
$ RCONDV, WORK, LWORK, RWORK, INFO )
+ implicit none
*
* -- LAPACK driver routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -312,8 +313,8 @@
LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE,
$ WNTSNN, WNTSNV
CHARACTER JOB, SIDE
- INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK,
- $ MINWRK, NOUT
+ INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K,
+ $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
COMPLEX*16 TMP
* ..
@@ -323,7 +324,7 @@
* ..
* .. External Subroutines ..
EXTERNAL DLABAD, DLASCL, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL,
- $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC,
+ $ ZGEHRD, ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC3,
$ ZTRSNA, ZUNGHR
* ..
* .. External Functions ..
@@ -387,9 +388,19 @@
MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
*
IF( WANTVL ) THEN
+ CALL ZTREVC3( 'L', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, LWORK_TREVC )
CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
$ WORK, -1, INFO )
ELSE IF( WANTVR ) THEN
+ CALL ZTREVC3( 'R', 'B', SELECT, N, A, LDA,
+ $ VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK, -1, RWORK, -1, IERR )
+ LWORK_TREVC = INT( WORK(1) )
+ MAXWRK = MAX( MAXWRK, LWORK_TREVC )
CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
$ WORK, -1, INFO )
ELSE
@@ -401,7 +412,7 @@
$ WORK, -1, INFO )
END IF
END IF
- HSWORK = WORK( 1 )
+ HSWORK = INT( WORK(1) )
*
IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
MINWRK = 2*N
@@ -567,11 +578,12 @@
IF( WANTVL .OR. WANTVR ) THEN
*
* Compute left and/or right eigenvectors
-* (CWorkspace: need 2*N)
+* (CWorkspace: need 2*N, prefer N + 2*N*NB)
* (RWorkspace: need N)
*
- CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
- $ N, NOUT, WORK( IWRK ), RWORK, IERR )
+ CALL ZTREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
+ $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1,
+ $ RWORK, N, IERR )
END IF
*
* Compute condition numbers if desired
diff --git a/SRC/ztrevc3.f b/SRC/ztrevc3.f
new file mode 100644
index 00000000..22654856
--- /dev/null
+++ b/SRC/ztrevc3.f
@@ -0,0 +1,630 @@
+*> \brief \b ZTREVC3
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZTREVC3 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrevc3.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrevc3.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrevc3.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
+* VR, LDVR, MM, M, WORK, LWORK, RWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER HOWMNY, SIDE
+* INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+* LOGICAL SELECT( * )
+* DOUBLE PRECISION RWORK( * )
+* COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+* $ WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> 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.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*> SIDE is CHARACTER*1
+*> = 'R': compute right eigenvectors only;
+*> = 'L': compute left eigenvectors only;
+*> = 'B': compute both right and left eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] HOWMNY
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[in] SELECT
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix T. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] T
+*> \verbatim
+*> T is COMPLEX*16 array, dimension (LDT,N)
+*> The upper triangular matrix T. T is modified, but restored
+*> on exit.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*> LDT is INTEGER
+*> The leading dimension of the array T. LDT >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] VL
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the array VL.
+*> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[in,out] VR
+*> \verbatim
+*> 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'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the array VR.
+*> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[in] MM
+*> \verbatim
+*> MM is INTEGER
+*> The number of columns in the arrays VL and/or VR. MM >= M.
+*> \endverbatim
+*>
+*> \param[out] M
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension (LRWORK)
+*> \endverbatim
+*>
+*> \param[in] LRWORK
+*> \verbatim
+*> 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.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+* @precisions fortran z -> c
+*
+*> \ingroup complex16OTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> 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|.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE ZTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
+ $ LDVR, MM, M, WORK, LWORK, RWORK, LRWORK, INFO)
+ IMPLICIT NONE
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER HOWMNY, SIDE
+ INTEGER INFO, LDT, LDVL, LDVR, LWORK, LRWORK, M, MM, N
+* ..
+* .. Array Arguments ..
+ LOGICAL SELECT( * )
+ DOUBLE PRECISION RWORK( * )
+ COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
+ $ WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ COMPLEX*16 CZERO, CONE
+ PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
+ $ CONE = ( 1.0D+0, 0.0D+0 ) )
+ INTEGER NBMIN, NBMAX
+ PARAMETER ( NBMIN = 8, NBMAX = 128 )
+* ..
+* .. Local Scalars ..
+ LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, RIGHTV, SOMEV
+ INTEGER I, II, IS, J, K, KI, IV, MAXWRK, NB
+ DOUBLE PRECISION OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
+ COMPLEX*16 CDUM
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ILAENV, IZAMAX
+ DOUBLE PRECISION DLAMCH, DZASUM
+ EXTERNAL LSAME, ILAENV, IZAMAX, DLAMCH, DZASUM
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, DCMPLX, CONJG, AIMAG, MAX
+* ..
+* .. Statement Functions ..
+ DOUBLE PRECISION CABS1
+* ..
+* .. Statement Function definitions ..
+ CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( AIMAG( CDUM ) )
+* ..
+* .. Executable Statements ..
+*
+* Decode and test the input parameters
+*
+ BOTHV = LSAME( SIDE, 'B' )
+ RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
+ LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
+*
+ ALLV = LSAME( HOWMNY, 'A' )
+ OVER = LSAME( HOWMNY, 'B' )
+ SOMEV = LSAME( HOWMNY, 'S' )
+*
+* Set M to the number of columns required to store the selected
+* eigenvectors.
+*
+ IF( SOMEV ) THEN
+ M = 0
+ DO 10 J = 1, N
+ IF( SELECT( J ) )
+ $ M = M + 1
+ 10 CONTINUE
+ ELSE
+ M = N
+ END IF
+*
+ INFO = 0
+ NB = ILAENV( 1, 'ZTREVC', SIDE // HOWMNY, N, -1, -1, -1 )
+ MAXWRK = N + 2*N*NB
+ WORK(1) = MAXWRK
+ RWORK(1) = N
+ LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
+ IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
+ INFO = -1
+ ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
+ INFO = -8
+ ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
+ INFO = -10
+ ELSE IF( MM.LT.M ) THEN
+ INFO = -11
+ ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -14
+ ELSE IF ( LRWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -16
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZTREVC3', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Use blocked version of back-transformation if sufficient workspace.
+* Zero-out the workspace to avoid potential NaN propagation.
+*
+ IF( OVER .AND. LWORK .GE. N + 2*N*NBMIN ) THEN
+ NB = (LWORK - N) / (2*N)
+ NB = MIN( NB, NBMAX )
+ CALL ZLASET( 'F', N, 1+2*NB, CZERO, CZERO, WORK, N )
+ ELSE
+ NB = 1
+ END IF
+*
+* Set the constants to control overflow.
+*
+ UNFL = DLAMCH( 'Safe minimum' )
+ OVFL = ONE / UNFL
+ CALL DLABAD( UNFL, OVFL )
+ ULP = DLAMCH( 'Precision' )
+ SMLNUM = UNFL*( N / ULP )
+*
+* Store the diagonal elements of T in working array WORK.
+*
+ DO 20 I = 1, N
+ WORK( I ) = T( I, I )
+ 20 CONTINUE
+*
+* Compute 1-norm of each column of strictly upper triangular
+* part of T to control overflow in triangular solver.
+*
+ RWORK( 1 ) = ZERO
+ DO 30 J = 2, N
+ RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
+ 30 CONTINUE
+*
+ IF( RIGHTV ) THEN
+*
+* ============================================================
+* Compute right eigenvectors.
+*
+* IV is index of column in current block.
+* Non-blocked version always uses IV=NB=1;
+* blocked version starts with IV=NB, goes down to 1.
+* (Note the "0-th" column is used to store the original diagonal.)
+ IV = NB
+ IS = M
+ DO 80 KI = N, 1, -1
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 80
+ END IF
+ SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
+*
+* --------------------------------------------------------
+* Complex right eigenvector
+*
+ WORK( KI + IV*N ) = CONE
+*
+* Form right-hand side.
+*
+ DO 40 K = 1, KI - 1
+ WORK( K + IV*N ) = -T( K, KI )
+ 40 CONTINUE
+*
+* Solve upper triangular system:
+* [ T(1:KI-1,1:KI-1) - T(KI,KI) ]*X = SCALE*WORK.
+*
+ DO 50 K = 1, KI - 1
+ T( K, K ) = T( K, K ) - T( KI, KI )
+ IF( CABS1( T( K, K ) ).LT.SMIN )
+ $ T( K, K ) = SMIN
+ 50 CONTINUE
+*
+ IF( KI.GT.1 ) THEN
+ CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
+ $ KI-1, T, LDT, WORK( 1 + IV*N ), SCALE,
+ $ RWORK, INFO )
+ WORK( KI + IV*N ) = SCALE
+ END IF
+*
+* Copy the vector x or Q*x to VR and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VR and normalize.
+ CALL ZCOPY( KI, WORK( 1 + IV*N ), 1, VR( 1, IS ), 1 )
+*
+ II = IZAMAX( KI, VR( 1, IS ), 1 )
+ REMAX = ONE / CABS1( VR( II, IS ) )
+ CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
+*
+ DO 60 K = KI + 1, N
+ VR( K, IS ) = CZERO
+ 60 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.GT.1 )
+ $ CALL ZGEMV( 'N', N, KI-1, CONE, VR, LDVR,
+ $ WORK( 1 + IV*N ), 1, DCMPLX( SCALE ),
+ $ VR( 1, KI ), 1 )
+*
+ II = IZAMAX( N, VR( 1, KI ), 1 )
+ REMAX = ONE / CABS1( VR( II, KI ) )
+ CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out below vector
+ DO K = KI + 1, N
+ WORK( K + IV*N ) = CZERO
+ END DO
+*
+* Columns IV:NB of work are valid vectors.
+* When the number of vectors stored reaches NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.EQ.1) .OR. (KI.EQ.1) ) THEN
+ CALL ZGEMM( 'N', 'N', N, NB-IV+1, KI+NB-IV, CONE,
+ $ VR, LDVR,
+ $ WORK( 1 + (IV)*N ), N,
+ $ CZERO,
+ $ WORK( 1 + (NB+IV)*N ), N )
+* normalize vectors
+ DO K = IV, NB
+ II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
+ CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL ZLACPY( 'F', N, NB-IV+1,
+ $ WORK( 1 + (NB+IV)*N ), N,
+ $ VR( 1, KI ), LDVR )
+ IV = NB
+ ELSE
+ IV = IV - 1
+ END IF
+ END IF
+*
+* Restore the original diagonal elements of T.
+*
+ DO 70 K = 1, KI - 1
+ T( K, K ) = WORK( K )
+ 70 CONTINUE
+*
+ IS = IS - 1
+ 80 CONTINUE
+ END IF
+*
+ IF( LEFTV ) THEN
+*
+* ============================================================
+* Compute left eigenvectors.
+*
+* IV is index of column in current block.
+* Non-blocked version always uses IV=1;
+* blocked version starts with IV=1, goes up to NB.
+* (Note the "0-th" column is used to store the original diagonal.)
+ IV = 1
+ IS = 1
+ DO 130 KI = 1, N
+*
+ IF( SOMEV ) THEN
+ IF( .NOT.SELECT( KI ) )
+ $ GO TO 130
+ END IF
+ SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
+*
+* --------------------------------------------------------
+* Complex left eigenvector
+*
+ WORK( KI + IV*N ) = CONE
+*
+* Form right-hand side.
+*
+ DO 90 K = KI + 1, N
+ WORK( K + IV*N ) = -CONJG( T( KI, K ) )
+ 90 CONTINUE
+*
+* Solve conjugate-transposed triangular system:
+* [ T(KI+1:N,KI+1:N) - T(KI,KI) ]**H * X = SCALE*WORK.
+*
+ DO 100 K = KI + 1, N
+ T( K, K ) = T( K, K ) - T( KI, KI )
+ IF( CABS1( T( K, K ) ).LT.SMIN )
+ $ T( K, K ) = SMIN
+ 100 CONTINUE
+*
+ IF( KI.LT.N ) THEN
+ CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
+ $ 'Y', N-KI, T( KI+1, KI+1 ), LDT,
+ $ WORK( KI+1 + IV*N ), SCALE, RWORK, INFO )
+ WORK( KI + IV*N ) = SCALE
+ END IF
+*
+* Copy the vector x or Q*x to VL and normalize.
+*
+ IF( .NOT.OVER ) THEN
+* ------------------------------
+* no back-transform: copy x to VL and normalize.
+ CALL ZCOPY( N-KI+1, WORK( KI + IV*N ), 1, VL(KI,IS), 1 )
+*
+ II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
+ REMAX = ONE / CABS1( VL( II, IS ) )
+ CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
+*
+ DO 110 K = 1, KI - 1
+ VL( K, IS ) = CZERO
+ 110 CONTINUE
+*
+ ELSE IF( NB.EQ.1 ) THEN
+* ------------------------------
+* version 1: back-transform each vector with GEMV, Q*x.
+ IF( KI.LT.N )
+ $ CALL ZGEMV( 'N', N, N-KI, CONE, VL( 1, KI+1 ), LDVL,
+ $ WORK( KI+1 + IV*N ), 1, DCMPLX( SCALE ),
+ $ VL( 1, KI ), 1 )
+*
+ II = IZAMAX( N, VL( 1, KI ), 1 )
+ REMAX = ONE / CABS1( VL( II, KI ) )
+ CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
+*
+ ELSE
+* ------------------------------
+* version 2: back-transform block of vectors with GEMM
+* zero out above vector
+* could go from KI-NV+1 to KI-1
+ DO K = 1, KI - 1
+ WORK( K + IV*N ) = CZERO
+ END DO
+*
+* Columns 1:IV of work are valid vectors.
+* When the number of vectors stored reaches NB,
+* or if this was last vector, do the GEMM
+ IF( (IV.EQ.NB) .OR. (KI.EQ.N) ) THEN
+ CALL ZGEMM( 'N', 'N', N, IV, N-KI+IV, ONE,
+ $ VL( 1, KI-IV+1 ), LDVL,
+ $ WORK( KI-IV+1 + (1)*N ), N,
+ $ CZERO,
+ $ WORK( 1 + (NB+1)*N ), N )
+* normalize vectors
+ DO K = 1, IV
+ II = IZAMAX( N, WORK( 1 + (NB+K)*N ), 1 )
+ REMAX = ONE / CABS1( WORK( II + (NB+K)*N ) )
+ CALL ZDSCAL( N, REMAX, WORK( 1 + (NB+K)*N ), 1 )
+ END DO
+ CALL ZLACPY( 'F', N, IV,
+ $ WORK( 1 + (NB+1)*N ), N,
+ $ VL( 1, KI-IV+1 ), LDVL )
+ IV = 1
+ ELSE
+ IV = IV + 1
+ END IF
+ END IF
+*
+* Restore the original diagonal elements of T.
+*
+ DO 120 K = KI + 1, N
+ T( K, K ) = WORK( K )
+ 120 CONTINUE
+*
+ IS = IS + 1
+ 130 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of ZTREVC3
+*
+ END
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-23 01:53:15 +0000