--- rpl/lapack/lapack/dsygvd.f 2010/08/07 13:22:27 1.5 +++ rpl/lapack/lapack/dsygvd.f 2023/08/07 08:39:08 1.19 @@ -1,10 +1,233 @@ +*> \brief \b DSYGVD +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DSYGVD + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, +* LWORK, IWORK, LIWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER JOBZ, UPLO +* INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N +* .. +* .. Array Arguments .. +* INTEGER IWORK( * ) +* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DSYGVD computes all the eigenvalues, and optionally, the eigenvectors +*> of a real generalized symmetric-definite eigenproblem, of the form +*> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and +*> B are assumed to be symmetric and B is also positive definite. +*> If eigenvectors are desired, it uses a divide and conquer algorithm. +*> +*> The divide and conquer algorithm makes very mild assumptions about +*> floating point arithmetic. It will work on machines with a guard +*> digit in add/subtract, or on those binary machines without guard +*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or +*> Cray-2. It could conceivably fail on hexadecimal or decimal machines +*> without guard digits, but we know of none. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] ITYPE +*> \verbatim +*> ITYPE is INTEGER +*> Specifies the problem type to be solved: +*> = 1: A*x = (lambda)*B*x +*> = 2: A*B*x = (lambda)*x +*> = 3: B*A*x = (lambda)*x +*> \endverbatim +*> +*> \param[in] JOBZ +*> \verbatim +*> JOBZ is CHARACTER*1 +*> = 'N': Compute eigenvalues only; +*> = 'V': Compute eigenvalues and eigenvectors. +*> \endverbatim +*> +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangles of A and B are stored; +*> = 'L': Lower triangles of A and B are stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrices A and B. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA, N) +*> On entry, the symmetric matrix A. If UPLO = 'U', the +*> leading N-by-N upper triangular part of A contains the +*> upper triangular part of the matrix A. If UPLO = 'L', +*> the leading N-by-N lower triangular part of A contains +*> the lower triangular part of the matrix A. +*> +*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the +*> matrix Z of eigenvectors. The eigenvectors are normalized +*> as follows: +*> if ITYPE = 1 or 2, Z**T*B*Z = I; +*> if ITYPE = 3, Z**T*inv(B)*Z = I. +*> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') +*> or the lower triangle (if UPLO='L') of A, including the +*> diagonal, is destroyed. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[in,out] B +*> \verbatim +*> B is DOUBLE PRECISION array, dimension (LDB, N) +*> On entry, the symmetric matrix B. If UPLO = 'U', the +*> leading N-by-N upper triangular part of B contains the +*> upper triangular part of the matrix B. If UPLO = 'L', +*> the leading N-by-N lower triangular part of B contains +*> the lower triangular part of the matrix B. +*> +*> On exit, if INFO <= N, the part of B containing the matrix is +*> overwritten by the triangular factor U or L from the Cholesky +*> factorization B = U**T*U or B = L*L**T. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is INTEGER +*> The leading dimension of the array B. LDB >= max(1,N). +*> \endverbatim +*> +*> \param[out] W +*> \verbatim +*> W is DOUBLE PRECISION array, dimension (N) +*> If INFO = 0, the eigenvalues in ascending order. +*> \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. +*> If N <= 1, LWORK >= 1. +*> If JOBZ = 'N' and N > 1, LWORK >= 2*N+1. +*> If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2. +*> +*> If LWORK = -1, then a workspace query is assumed; the routine +*> only calculates the optimal sizes of the WORK and IWORK +*> arrays, returns these values as the first entries of the WORK +*> and IWORK arrays, and no error message related to LWORK or +*> LIWORK is issued by XERBLA. +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) +*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. +*> \endverbatim +*> +*> \param[in] LIWORK +*> \verbatim +*> LIWORK is INTEGER +*> The dimension of the array IWORK. +*> If N <= 1, LIWORK >= 1. +*> If JOBZ = 'N' and N > 1, LIWORK >= 1. +*> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. +*> +*> If LIWORK = -1, then a workspace query is assumed; the +*> routine only calculates the optimal sizes of the WORK and +*> IWORK arrays, returns these values as the first entries of +*> the WORK and IWORK arrays, and no error message related to +*> LWORK or LIWORK 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 +*> > 0: DPOTRF or DSYEVD returned an error code: +*> <= N: if INFO = i and JOBZ = 'N', then the algorithm +*> failed to converge; i off-diagonal elements of an +*> intermediate tridiagonal form did not converge to +*> zero; +*> if INFO = i and JOBZ = 'V', then the algorithm +*> failed to compute an eigenvalue while working on +*> the submatrix lying in rows and columns INFO/(N+1) +*> through mod(INFO,N+1); +*> > N: if INFO = N + i, for 1 <= i <= N, then the leading +*> minor of order i of B is not positive definite. +*> The factorization of B could not be completed and +*> no eigenvalues or eigenvectors were computed. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup doubleSYeigen +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> Modified so that no backsubstitution is performed if DSYEVD fails to +*> converge (NEIG in old code could be greater than N causing out of +*> bounds reference to A - reported by Ralf Meyer). Also corrected the +*> description of INFO and the test on ITYPE. Sven, 16 Feb 05. +*> \endverbatim +* +*> \par Contributors: +* ================== +*> +*> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA +*> +* ===================================================================== SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, $ LWORK, IWORK, LIWORK, 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 JOBZ, UPLO @@ -15,135 +238,6 @@ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) * .. * -* Purpose -* ======= -* -* DSYGVD computes all the eigenvalues, and optionally, the eigenvectors -* of a real generalized symmetric-definite eigenproblem, of the form -* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and -* B are assumed to be symmetric and B is also positive definite. -* If eigenvectors are desired, it uses a divide and conquer algorithm. -* -* The divide and conquer algorithm makes very mild assumptions about -* floating point arithmetic. It will work on machines with a guard -* digit in add/subtract, or on those binary machines without guard -* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or -* Cray-2. It could conceivably fail on hexadecimal or decimal machines -* without guard digits, but we know of none. -* -* Arguments -* ========= -* -* ITYPE (input) INTEGER -* Specifies the problem type to be solved: -* = 1: A*x = (lambda)*B*x -* = 2: A*B*x = (lambda)*x -* = 3: B*A*x = (lambda)*x -* -* JOBZ (input) CHARACTER*1 -* = 'N': Compute eigenvalues only; -* = 'V': Compute eigenvalues and eigenvectors. -* -* UPLO (input) CHARACTER*1 -* = 'U': Upper triangles of A and B are stored; -* = 'L': Lower triangles of A and B are stored. -* -* N (input) INTEGER -* The order of the matrices A and B. N >= 0. -* -* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) -* On entry, the symmetric matrix A. If UPLO = 'U', the -* leading N-by-N upper triangular part of A contains the -* upper triangular part of the matrix A. If UPLO = 'L', -* the leading N-by-N lower triangular part of A contains -* the lower triangular part of the matrix A. -* -* On exit, if JOBZ = 'V', then if INFO = 0, A contains the -* matrix Z of eigenvectors. The eigenvectors are normalized -* as follows: -* if ITYPE = 1 or 2, Z**T*B*Z = I; -* if ITYPE = 3, Z**T*inv(B)*Z = I. -* If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') -* or the lower triangle (if UPLO='L') of A, including the -* diagonal, is destroyed. -* -* LDA (input) INTEGER -* The leading dimension of the array A. LDA >= max(1,N). -* -* B (input/output) DOUBLE PRECISION array, dimension (LDB, N) -* On entry, the symmetric matrix B. If UPLO = 'U', the -* leading N-by-N upper triangular part of B contains the -* upper triangular part of the matrix B. If UPLO = 'L', -* the leading N-by-N lower triangular part of B contains -* the lower triangular part of the matrix B. -* -* On exit, if INFO <= N, the part of B containing the matrix is -* overwritten by the triangular factor U or L from the Cholesky -* factorization B = U**T*U or B = L*L**T. -* -* LDB (input) INTEGER -* The leading dimension of the array B. LDB >= max(1,N). -* -* W (output) DOUBLE PRECISION array, dimension (N) -* If INFO = 0, the eigenvalues in ascending order. -* -* 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. -* If N <= 1, LWORK >= 1. -* If JOBZ = 'N' and N > 1, LWORK >= 2*N+1. -* If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2. -* -* If LWORK = -1, then a workspace query is assumed; the routine -* only calculates the optimal sizes of the WORK and IWORK -* arrays, returns these values as the first entries of the WORK -* and IWORK arrays, and no error message related to LWORK or -* LIWORK is issued by XERBLA. -* -* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) -* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. -* -* LIWORK (input) INTEGER -* The dimension of the array IWORK. -* If N <= 1, LIWORK >= 1. -* If JOBZ = 'N' and N > 1, LIWORK >= 1. -* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. -* -* If LIWORK = -1, then a workspace query is assumed; the -* routine only calculates the optimal sizes of the WORK and -* IWORK arrays, returns these values as the first entries of -* the WORK and IWORK arrays, and no error message related to -* LWORK or LIWORK is issued by XERBLA. -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* > 0: DPOTRF or DSYEVD returned an error code: -* <= N: if INFO = i and JOBZ = 'N', then the algorithm -* failed to converge; i off-diagonal elements of an -* intermediate tridiagonal form did not converge to -* zero; -* if INFO = i and JOBZ = 'V', then the algorithm -* failed to compute an eigenvalue while working on -* the submatrix lying in rows and columns INFO/(N+1) -* through mod(INFO,N+1); -* > N: if INFO = N + i, for 1 <= i <= N, then the leading -* minor of order i of B is not positive definite. -* The factorization of B could not be completed and -* no eigenvalues or eigenvectors were computed. -* -* Further Details -* =============== -* -* Based on contributions by -* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA -* -* Modified so that no backsubstitution is performed if DSYEVD fails to -* converge (NEIG in old code could be greater than N causing out of -* bounds reference to A - reported by Ralf Meyer). Also corrected the -* description of INFO and the test on ITYPE. Sven, 16 Feb 05. * ===================================================================== * * .. Parameters .. @@ -236,8 +330,8 @@ CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK, $ INFO ) - LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) - LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) + LOPT = INT( MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) ) + LIOPT = INT( MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) ) * IF( WANTZ .AND. INFO.EQ.0 ) THEN * @@ -246,7 +340,7 @@ IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN * * For A*x=(lambda)*B*x and A*B*x=(lambda)*x; -* backtransform eigenvectors: x = inv(L)'*y or inv(U)*y +* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y * IF( UPPER ) THEN TRANS = 'N' @@ -260,7 +354,7 @@ ELSE IF( ITYPE.EQ.3 ) THEN * * For B*A*x=(lambda)*x; -* backtransform eigenvectors: x = L*y or U'*y +* backtransform eigenvectors: x = L*y or U**T*y * IF( UPPER ) THEN TRANS = 'T'