From ed2ea1af894955ddd1ddfd0acb15e1c07d459f1e Mon Sep 17 00:00:00 2001 From: Julie Date: Sun, 12 Jun 2016 23:21:04 -0700 Subject: 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.) --- SRC/CMakeLists.txt | 8 +- SRC/Makefile | 8 +- SRC/cgeev.f | 31 +- SRC/cgeevx.f | 26 +- SRC/ctrevc3.f | 630 +++++++++++++++++++++++++ SRC/dgeev.f | 33 +- SRC/dgeevx.f | 29 +- SRC/dtrevc3.f | 1303 ++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/ilaenv.f | 6 + SRC/sgeev.f | 33 +- SRC/sgeevx.f | 30 +- SRC/strevc3.f | 1303 ++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/zgeev.f | 30 +- SRC/zgeevx.f | 26 +- SRC/ztrevc3.f | 630 +++++++++++++++++++++++++ 15 files changed, 4045 insertions(+), 81 deletions(-) create mode 100644 SRC/ctrevc3.f create mode 100644 SRC/dtrevc3.f create mode 100644 SRC/strevc3.f create mode 100644 SRC/ztrevc3.f 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 +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \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 -- cgit v1.2.3