dgges.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 550 行 · 第 1/2 页
C
550 行
#include "f2c.h"
#include "netlib.h"
extern double sqrt(double); /* #include <math.h> */
/* Table of constant values */
static integer c__1 = 1;
static integer c__0 = 0;
static integer c_n1 = -1;
static doublereal c_b33 = 0.;
static doublereal c_b34 = 1.;
/* Subroutine */ void dgges_(jobvsl, jobvsr, sort, delctg, n, a, lda, b, ldb,
sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
char *jobvsl, *jobvsr, *sort;
logical (*delctg) (doublereal*,doublereal*,doublereal*);
integer *n;
doublereal *a;
integer *lda;
doublereal *b;
integer *ldb, *sdim;
doublereal *alphar, *alphai, *beta, *vsl;
integer *ldvsl;
doublereal *vsr;
integer *ldvsr;
doublereal *work;
integer *lwork;
logical *bwork;
integer *info;
{
/* System generated locals */
integer i__1;
/* Local variables */
static doublereal anrm, bnrm;
static integer idum[1], ierr, itau, iwrk;
static doublereal pvsl, pvsr;
static integer i;
static integer ileft, icols;
static logical cursl, ilvsl, ilvsr;
static integer irows;
static logical lst2sl;
static integer ip;
static logical ilascl, ilbscl;
static doublereal safmin;
static doublereal safmax;
static doublereal bignum;
static integer ijobvl, iright, ijobvr;
static doublereal anrmto, bnrmto;
static logical lastsl;
static integer minwrk, maxwrk;
static doublereal smlnum;
static logical wantst, lquery;
static doublereal dif[2];
static integer ihi, ilo;
static doublereal eps;
/* -- LAPACK driver routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* June 30, 1999 */
/* Purpose */
/* ======= */
/* */
/* DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B), */
/* the generalized eigenvalues, the generalized real Schur form (S,T), */
/* optionally, the left and/or right matrices of Schur vectors (VSL and */
/* VSR). This gives the generalized Schur factorization */
/* */
/* (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T ) */
/* */
/* Optionally, it also orders the eigenvalues so that a selected cluster */
/* of eigenvalues appears in the leading diagonal blocks of the upper */
/* quasi-triangular matrix S and the upper triangular matrix T.The */
/* leading columns of VSL and VSR then form an orthonormal basis for the */
/* corresponding left and right eigenspaces (deflating subspaces). */
/* */
/* (If only the generalized eigenvalues are needed, use the driver */
/* DGGEV instead, which is faster.) */
/* */
/* A generalized eigenvalue for a pair of matrices (A,B) is a scalar w */
/* or a ratio alpha/beta = w, such that A - w*B is singular. It is */
/* usually represented as the pair (alpha,beta), as there is a */
/* reasonable interpretation for beta=0 or both being zero. */
/* */
/* A pair of matrices (S,T) is in generalized real Schur form if T is */
/* upper triangular with non-negative diagonal and S is block upper */
/* triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond */
/* to real generalized eigenvalues, while 2-by-2 blocks of S will be */
/* "standardized" by making the corresponding elements of T have the */
/* form: */
/* [ a 0 ] */
/* [ 0 b ] */
/* */
/* and the pair of corresponding 2-by-2 blocks in S and T will have a */
/* complex conjugate pair of generalized eigenvalues. */
/* */
/* */
/* Arguments */
/* ========= */
/* */
/* JOBVSL (input) CHARACTER*1 */
/* = 'N': do not compute the left Schur vectors; */
/* = 'V': compute the left Schur vectors. */
/* */
/* JOBVSR (input) CHARACTER*1 */
/* = 'N': do not compute the right Schur vectors; */
/* = 'V': compute the right Schur vectors. */
/* */
/* SORT (input) CHARACTER*1 */
/* Specifies whether or not to order the eigenvalues on the */
/* diagonal of the generalized Schur form. */
/* = 'N': Eigenvalues are not ordered; */
/* = 'S': Eigenvalues are ordered (see DELZTG); */
/* */
/* DELZTG (input) LOGICAL FUNCTION of three DOUBLE PRECISION arguments */
/* DELZTG must be declared EXTERNAL in the calling subroutine. */
/* If SORT = 'N', DELZTG is not referenced. */
/* If SORT = 'S', DELZTG is used to select eigenvalues to sort */
/* to the top left of the Schur form. */
/* An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if */
/* DELZTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either */
/* one of a complex conjugate pair of eigenvalues is selected, */
/* then both complex eigenvalues are selected. */
/* */
/* Note that in the ill-conditioned case, a selected complex */
/* eigenvalue may no longer satisfy DELZTG(ALPHAR(j),ALPHAI(j), */
/* BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2 */
/* in this case. */
/* */
/* N (input) INTEGER */
/* The order of the matrices A, B, VSL, and VSR. N >= 0. */
/* */
/* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
/* On entry, the first of the pair of matrices. */
/* On exit, A has been overwritten by its generalized Schur */
/* form S. */
/* */
/* LDA (input) INTEGER */
/* The leading dimension of A. LDA >= max(1,N). */
/* */
/* B (input/output) DOUBLE PRECISION array, dimension (LDB, N) */
/* On entry, the second of the pair of matrices. */
/* On exit, B has been overwritten by its generalized Schur */
/* form T. */
/* */
/* LDB (input) INTEGER */
/* The leading dimension of B. LDB >= max(1,N). */
/* */
/* SDIM (output) INTEGER */
/* If SORT = 'N', SDIM = 0. */
/* If SORT = 'S', SDIM = number of eigenvalues (after sorting) */
/* for which DELZTG is true. (Complex conjugate pairs for which */
/* DELZTG is true for either eigenvalue count as 2.) */
/* */
/* ALPHAR (output) DOUBLE PRECISION array, dimension (N) */
/* ALPHAI (output) DOUBLE PRECISION array, dimension (N) */
/* BETA (output) DOUBLE PRECISION array, dimension (N) */
/* On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will */
/* be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i, */
/* and BETA(j),j=1,...,N are the diagonals of the complex Schur */
/* form (S,T) that would result if the 2-by-2 diagonal blocks of */
/* the real Schur form of (A,B) were further reduced to */
/* triangular form using 2-by-2 complex unitary transformations. */
/* 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) negative. */
/* */
/* Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j) */
/* may easily over- or underflow, and BETA(j) may even be zero. */
/* Thus, the user should avoid naively computing the ratio. */
/* However, ALPHAR and ALPHAI will be always less than and */
/* usually comparable with norm(A) in magnitude, and BETA always */
/* less than and usually comparable with norm(B). */
/* */
/* VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N) */
/* If JOBVSL = 'V', VSL will contain the left Schur vectors. */
/* Not referenced if JOBVSL = 'N'. */
/* */
/* LDVSL (input) INTEGER */
/* The leading dimension of the matrix VSL. LDVSL >=1, and */
/* if JOBVSL = 'V', LDVSL >= N. */
/* */
/* VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N) */
/* If JOBVSR = 'V', VSR will contain the right Schur vectors. */
/* Not referenced if JOBVSR = 'N'. */
/* */
/* LDVSR (input) INTEGER */
/* The leading dimension of the matrix VSR. LDVSR >= 1, and */
/* if JOBVSR = 'V', LDVSR >= N. */
/* */
/* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
/* */
/* LWORK (input) INTEGER */
/* The dimension of the array WORK. LWORK >= 8*N+16. */
/* */
/* 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. */
/* */
/* BWORK (workspace) LOGICAL array, dimension (N) */
/* Not referenced if SORT = 'N'. */
/* */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value. */
/* = 1,...,N: */
/* The QZ iteration failed. (A,B) are not in Schur */
/* form, but ALPHAR(j), ALPHAI(j), and BETA(j) should */
/* be correct for j=INFO+1,...,N. */
/* > N: =N+1: other than QZ iteration failed in DHGEQZ. */
/* =N+2: after reordering, roundoff changed values of */
/* some complex eigenvalues so that leading */
/* eigenvalues in the Generalized Schur form no */
/* longer satisfy DELZTG=.TRUE. This could also */
/* be caused due to scaling. */
/* =N+3: reordering failed in DTGSEN. */
/* */
/* ===================================================================== */
/* Decode the input arguments */
if (lsame_(jobvsl, "N")) {
ijobvl = 1;
ilvsl = FALSE_;
} else if (lsame_(jobvsl, "V")) {
ijobvl = 2;
ilvsl = TRUE_;
} else {
ijobvl = -1;
ilvsl = FALSE_;
}
if (lsame_(jobvsr, "N")) {
ijobvr = 1;
ilvsr = FALSE_;
} else if (lsame_(jobvsr, "V")) {
ijobvr = 2;
ilvsr = TRUE_;
} else {
ijobvr = -1;
ilvsr = FALSE_;
}
wantst = lsame_(sort, "S");
/* Test the input arguments */
*info = 0;
lquery = *lwork == -1;
if (ijobvl <= 0) {
*info = -1;
} else if (ijobvr <= 0) {
*info = -2;
} else if (! wantst && ! lsame_(sort, "N")) {
*info = -3;
} else if (*n < 0) {
*info = -5;
} else if (*lda < max(1,*n)) {
*info = -7;
} else if (*ldb < max(1,*n)) {
*info = -9;
} else if (*ldvsl < 1 || (ilvsl && *ldvsl < *n)) {
*info = -15;
} else if (*ldvsr < 1 || (ilvsr && *ldvsr < *n)) {
*info = -17;
}
/* Compute workspace */
/* (Note: Comments in the code beginning "Workspace:" describe the */
/* minimal amount of workspace needed at that point in the code, */
/* as well as the preferred amount for good performance. */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?