--- rpl/lapack/lapack/dgegv.f 2010/08/07 13:22:13 1.5
+++ rpl/lapack/lapack/dgegv.f 2023/08/07 08:38:48 1.17
@@ -1,10 +1,312 @@
+*> \brief DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGEGV + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
+* BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER JOBVL, JOBVR
+* INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
+* $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
+* $ VR( LDVR, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> This routine is deprecated and has been replaced by routine DGGEV.
+*>
+*> DGEGV computes the eigenvalues and, optionally, the left and/or right
+*> eigenvectors of a real matrix pair (A,B).
+*> Given two square matrices A and B,
+*> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
+*> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
+*> that
+*>
+*> A*x = lambda*B*x.
+*>
+*> An alternate form is to find the eigenvalues mu and corresponding
+*> eigenvectors y such that
+*>
+*> mu*A*y = B*y.
+*>
+*> These two forms are equivalent with mu = 1/lambda and x = y if
+*> neither lambda nor mu is zero. In order to deal with the case that
+*> lambda or mu is zero or small, two values alpha and beta are returned
+*> for each eigenvalue, such that lambda = alpha/beta and
+*> mu = beta/alpha.
+*>
+*> The vectors x and y in the above equations are right eigenvectors of
+*> the matrix pair (A,B). Vectors u and v satisfying
+*>
+*> u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
+*>
+*> are left eigenvectors of (A,B).
+*>
+*> Note: this routine performs "full balancing" on A and B
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] JOBVL
+*> \verbatim
+*> JOBVL is CHARACTER*1
+*> = 'N': do not compute the left generalized eigenvectors;
+*> = 'V': compute the left generalized eigenvectors (returned
+*> in VL).
+*> \endverbatim
+*>
+*> \param[in] JOBVR
+*> \verbatim
+*> JOBVR is CHARACTER*1
+*> = 'N': do not compute the right generalized eigenvectors;
+*> = 'V': compute the right generalized eigenvectors (returned
+*> in VR).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A, B, VL, and VR. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array, dimension (LDA, N)
+*> On entry, the matrix A.
+*> If JOBVL = 'V' or JOBVR = 'V', then on exit A
+*> contains the real Schur form of A from the generalized Schur
+*> factorization of the pair (A,B) after balancing.
+*> If no eigenvectors were computed, then only the diagonal
+*> blocks from the Schur form will be correct. See DGGHRD and
+*> DHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is DOUBLE PRECISION array, dimension (LDB, N)
+*> On entry, the matrix B.
+*> If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
+*> upper triangular matrix obtained from B in the generalized
+*> Schur factorization of the pair (A,B) after balancing.
+*> If no eigenvectors were computed, then only those elements of
+*> B corresponding to the diagonal blocks from the Schur form of
+*> A will be correct. See DGGHRD and DHGEQZ for details.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] ALPHAR
+*> \verbatim
+*> ALPHAR is DOUBLE PRECISION array, dimension (N)
+*> The real parts of each scalar alpha defining an eigenvalue of
+*> GNEP.
+*> \endverbatim
+*>
+*> \param[out] ALPHAI
+*> \verbatim
+*> ALPHAI is DOUBLE PRECISION array, dimension (N)
+*> The imaginary parts of each scalar alpha defining an
+*> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
+*> eigenvalue is real; if positive, then the j-th and
+*> (j+1)-st eigenvalues are a complex conjugate pair, with
+*> ALPHAI(j+1) = -ALPHAI(j).
+*> \endverbatim
+*>
+*> \param[out] BETA
+*> \verbatim
+*> BETA is DOUBLE PRECISION array, dimension (N)
+*> The scalars beta that define the eigenvalues of GNEP.
+*>
+*> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
+*> beta = BETA(j) represent the j-th eigenvalue of the matrix
+*> pair (A,B), in one of the forms lambda = alpha/beta or
+*> mu = beta/alpha. Since either lambda or mu may overflow,
+*> they should not, in general, be computed.
+*> \endverbatim
+*>
+*> \param[out] VL
+*> \verbatim
+*> VL is DOUBLE PRECISION array, dimension (LDVL,N)
+*> If JOBVL = 'V', the left eigenvectors u(j) are stored
+*> in the columns of VL, in the same order as their eigenvalues.
+*> If the j-th eigenvalue is real, then u(j) = VL(:,j).
+*> If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*> pair, then
+*> u(j) = VL(:,j) + i*VL(:,j+1)
+*> and
+*> u(j+1) = VL(:,j) - i*VL(:,j+1).
+*>
+*> Each eigenvector is scaled so that its largest component has
+*> abs(real part) + abs(imag. part) = 1, except for eigenvectors
+*> corresponding to an eigenvalue with alpha = beta = 0, which
+*> are set to zero.
+*> Not referenced if JOBVL = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVL
+*> \verbatim
+*> LDVL is INTEGER
+*> The leading dimension of the matrix VL. LDVL >= 1, and
+*> if JOBVL = 'V', LDVL >= N.
+*> \endverbatim
+*>
+*> \param[out] VR
+*> \verbatim
+*> VR is DOUBLE PRECISION array, dimension (LDVR,N)
+*> If JOBVR = 'V', the right eigenvectors x(j) are stored
+*> in the columns of VR, in the same order as their eigenvalues.
+*> If the j-th eigenvalue is real, then x(j) = VR(:,j).
+*> If the j-th and (j+1)-st eigenvalues form a complex conjugate
+*> pair, then
+*> x(j) = VR(:,j) + i*VR(:,j+1)
+*> and
+*> x(j+1) = VR(:,j) - i*VR(:,j+1).
+*>
+*> Each eigenvector is scaled so that its largest component has
+*> abs(real part) + abs(imag. part) = 1, except for eigenvalues
+*> corresponding to an eigenvalue with alpha = beta = 0, which
+*> are set to zero.
+*> Not referenced if JOBVR = 'N'.
+*> \endverbatim
+*>
+*> \param[in] LDVR
+*> \verbatim
+*> LDVR is INTEGER
+*> The leading dimension of the matrix VR. LDVR >= 1, and
+*> if JOBVR = 'V', LDVR >= N.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK. LWORK >= max(1,8*N).
+*> For good performance, LWORK must generally be larger.
+*> To compute the optimal value of LWORK, call ILAENV to get
+*> blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
+*> NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
+*> The optimal LWORK is:
+*> 2*N + MAX( 6*N, N*(NB+1) ).
+*>
+*> 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.
+*> = 1,...,N:
+*> The QZ iteration failed. No eigenvectors have been
+*> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
+*> should be correct for j=INFO+1,...,N.
+*> > N: errors that usually indicate LAPACK problems:
+*> =N+1: error return from DGGBAL
+*> =N+2: error return from DGEQRF
+*> =N+3: error return from DORMQR
+*> =N+4: error return from DORGQR
+*> =N+5: error return from DGGHRD
+*> =N+6: error return from DHGEQZ (other than failed
+*> iteration)
+*> =N+7: error return from DTGEVC
+*> =N+8: error return from DGGBAK (computing VL)
+*> =N+9: error return from DGGBAK (computing VR)
+*> =N+10: error return from DLASCL (various calls)
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleGEeigen
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> Balancing
+*> ---------
+*>
+*> This driver calls DGGBAL to both permute and scale rows and columns
+*> of A and B. The permutations PL and PR are chosen so that PL*A*PR
+*> and PL*B*R will be upper triangular except for the diagonal blocks
+*> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
+*> possible. The diagonal scaling matrices DL and DR are chosen so
+*> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
+*> one (except for the elements that start out zero.)
+*>
+*> After the eigenvalues and eigenvectors of the balanced matrices
+*> have been computed, DGGBAK transforms the eigenvectors back to what
+*> they would have been (in perfect arithmetic) if they had not been
+*> balanced.
+*>
+*> Contents of A and B on Exit
+*> -------- -- - --- - -- ----
+*>
+*> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
+*> both), then on exit the arrays A and B will contain the real Schur
+*> form[*] of the "balanced" versions of A and B. If no eigenvectors
+*> are computed, then only the diagonal blocks will be correct.
+*>
+*> [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
+*> by Golub & van Loan, pub. by Johns Hopkins U. Press.
+*> \endverbatim
+*>
+* =====================================================================
SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
$ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
*
-* -- LAPACK driver routine (version 3.2) --
+* -- LAPACK driver routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-* November 2006
*
* .. Scalar Arguments ..
CHARACTER JOBVL, JOBVR
@@ -16,208 +318,6 @@
$ VR( LDVR, * ), WORK( * )
* ..
*
-* Purpose
-* =======
-*
-* This routine is deprecated and has been replaced by routine DGGEV.
-*
-* DGEGV computes the eigenvalues and, optionally, the left and/or right
-* eigenvectors of a real matrix pair (A,B).
-* Given two square matrices A and B,
-* the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
-* eigenvalues lambda and corresponding (non-zero) eigenvectors x such
-* that
-*
-* A*x = lambda*B*x.
-*
-* An alternate form is to find the eigenvalues mu and corresponding
-* eigenvectors y such that
-*
-* mu*A*y = B*y.
-*
-* These two forms are equivalent with mu = 1/lambda and x = y if
-* neither lambda nor mu is zero. In order to deal with the case that
-* lambda or mu is zero or small, two values alpha and beta are returned
-* for each eigenvalue, such that lambda = alpha/beta and
-* mu = beta/alpha.
-*
-* The vectors x and y in the above equations are right eigenvectors of
-* the matrix pair (A,B). Vectors u and v satisfying
-*
-* u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
-*
-* are left eigenvectors of (A,B).
-*
-* Note: this routine performs "full balancing" on A and B -- see
-* "Further Details", below.
-*
-* Arguments
-* =========
-*
-* JOBVL (input) CHARACTER*1
-* = 'N': do not compute the left generalized eigenvectors;
-* = 'V': compute the left generalized eigenvectors (returned
-* in VL).
-*
-* JOBVR (input) CHARACTER*1
-* = 'N': do not compute the right generalized eigenvectors;
-* = 'V': compute the right generalized eigenvectors (returned
-* in VR).
-*
-* N (input) INTEGER
-* The order of the matrices A, B, VL, and VR. N >= 0.
-*
-* A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
-* On entry, the matrix A.
-* If JOBVL = 'V' or JOBVR = 'V', then on exit A
-* contains the real Schur form of A from the generalized Schur
-* factorization of the pair (A,B) after balancing.
-* If no eigenvectors were computed, then only the diagonal
-* blocks from the Schur form will be correct. See DGGHRD and
-* DHGEQZ for details.
-*
-* LDA (input) INTEGER
-* The leading dimension of A. LDA >= max(1,N).
-*
-* B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
-* On entry, the matrix B.
-* If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
-* upper triangular matrix obtained from B in the generalized
-* Schur factorization of the pair (A,B) after balancing.
-* If no eigenvectors were computed, then only those elements of
-* B corresponding to the diagonal blocks from the Schur form of
-* A will be correct. See DGGHRD and DHGEQZ for details.
-*
-* LDB (input) INTEGER
-* The leading dimension of B. LDB >= max(1,N).
-*
-* ALPHAR (output) DOUBLE PRECISION array, dimension (N)
-* The real parts of each scalar alpha defining an eigenvalue of
-* GNEP.
-*
-* ALPHAI (output) DOUBLE PRECISION array, dimension (N)
-* The imaginary parts of each scalar alpha defining an
-* eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
-* eigenvalue is real; if positive, then the j-th and
-* (j+1)-st eigenvalues are a complex conjugate pair, with
-* ALPHAI(j+1) = -ALPHAI(j).
-*
-* BETA (output) DOUBLE PRECISION array, dimension (N)
-* The scalars beta that define the eigenvalues of GNEP.
-*
-* Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
-* beta = BETA(j) represent the j-th eigenvalue of the matrix
-* pair (A,B), in one of the forms lambda = alpha/beta or
-* mu = beta/alpha. Since either lambda or mu may overflow,
-* they should not, in general, be computed.
-*
-* VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
-* If JOBVL = 'V', the left eigenvectors u(j) are stored
-* in the columns of VL, in the same order as their eigenvalues.
-* If the j-th eigenvalue is real, then u(j) = VL(:,j).
-* If the j-th and (j+1)-st eigenvalues form a complex conjugate
-* pair, then
-* u(j) = VL(:,j) + i*VL(:,j+1)
-* and
-* u(j+1) = VL(:,j) - i*VL(:,j+1).
-*
-* Each eigenvector is scaled so that its largest component has
-* abs(real part) + abs(imag. part) = 1, except for eigenvectors
-* corresponding to an eigenvalue with alpha = beta = 0, which
-* are set to zero.
-* Not referenced if JOBVL = 'N'.
-*
-* LDVL (input) INTEGER
-* The leading dimension of the matrix VL. LDVL >= 1, and
-* if JOBVL = 'V', LDVL >= N.
-*
-* VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
-* If JOBVR = 'V', the right eigenvectors x(j) are stored
-* in the columns of VR, in the same order as their eigenvalues.
-* If the j-th eigenvalue is real, then x(j) = VR(:,j).
-* If the j-th and (j+1)-st eigenvalues form a complex conjugate
-* pair, then
-* x(j) = VR(:,j) + i*VR(:,j+1)
-* and
-* x(j+1) = VR(:,j) - i*VR(:,j+1).
-*
-* Each eigenvector is scaled so that its largest component has
-* abs(real part) + abs(imag. part) = 1, except for eigenvalues
-* corresponding to an eigenvalue with alpha = beta = 0, which
-* are set to zero.
-* Not referenced if JOBVR = 'N'.
-*
-* LDVR (input) INTEGER
-* The leading dimension of the matrix VR. LDVR >= 1, and
-* if JOBVR = 'V', LDVR >= N.
-*
-* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
-* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
-*
-* LWORK (input) INTEGER
-* The dimension of the array WORK. LWORK >= max(1,8*N).
-* For good performance, LWORK must generally be larger.
-* To compute the optimal value of LWORK, call ILAENV to get
-* blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
-* NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
-* The optimal LWORK is:
-* 2*N + MAX( 6*N, N*(NB+1) ).
-*
-* 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.
-*
-* INFO (output) INTEGER
-* = 0: successful exit
-* < 0: if INFO = -i, the i-th argument had an illegal value.
-* = 1,...,N:
-* The QZ iteration failed. No eigenvectors have been
-* calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
-* should be correct for j=INFO+1,...,N.
-* > N: errors that usually indicate LAPACK problems:
-* =N+1: error return from DGGBAL
-* =N+2: error return from DGEQRF
-* =N+3: error return from DORMQR
-* =N+4: error return from DORGQR
-* =N+5: error return from DGGHRD
-* =N+6: error return from DHGEQZ (other than failed
-* iteration)
-* =N+7: error return from DTGEVC
-* =N+8: error return from DGGBAK (computing VL)
-* =N+9: error return from DGGBAK (computing VR)
-* =N+10: error return from DLASCL (various calls)
-*
-* Further Details
-* ===============
-*
-* Balancing
-* ---------
-*
-* This driver calls DGGBAL to both permute and scale rows and columns
-* of A and B. The permutations PL and PR are chosen so that PL*A*PR
-* and PL*B*R will be upper triangular except for the diagonal blocks
-* A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
-* possible. The diagonal scaling matrices DL and DR are chosen so
-* that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
-* one (except for the elements that start out zero.)
-*
-* After the eigenvalues and eigenvectors of the balanced matrices
-* have been computed, DGGBAK transforms the eigenvectors back to what
-* they would have been (in perfect arithmetic) if they had not been
-* balanced.
-*
-* Contents of A and B on Exit
-* -------- -- - --- - -- ----
-*
-* If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
-* both), then on exit the arrays A and B will contain the real Schur
-* form[*] of the "balanced" versions of A and B. If no eigenvectors
-* are computed, then only the diagonal blocks will be correct.
-*
-* [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
-* by Golub & van Loan, pub. by Johns Hopkins U. Press.
-*
* =====================================================================
*
* .. Parameters ..