ztgevc.c
来自「算断裂的」· C语言 代码 · 共 948 行 · 第 1/2 页
C
948 行
#include "f2c.h"
/* Subroutine */ int ztgevc_(char *side, char *howmny, logical *select,
integer *n, doublecomplex *a, integer *lda, doublecomplex *b, integer
*ldb, doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *
ldvr, integer *mm, integer *m, doublecomplex *work, doublereal *rwork,
integer *info)
{
/* -- LAPACK routine (version 2.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
September 30, 1994
Purpose
=======
ZTGEVC computes some or all of the right and/or left generalized
eigenvectors of a pair of complex upper triangular matrices (A,B).
The right generalized eigenvector x and the left generalized
eigenvector y of (A,B) corresponding to a generalized eigenvalue
w are defined by:
(A - wB) * x = 0 and y**H * (A - wB) = 0
where y**H denotes the conjugate tranpose of y.
If an eigenvalue w is determined by zero diagonal elements of both A
and B, a unit vector is returned as the corresponding eigenvector.
If all eigenvectors are requested, the routine may either return
the matrices X and/or Y of right or left eigenvectors of (A,B), or
the products Z*X and/or Q*Y, where Z and Q are input unitary
matrices. If (A,B) was obtained from the generalized Schur
factorization of an original pair of matrices
(A0,B0) = (Q*A*Z**H,Q*B*Z**H),
then Z*X and Q*Y are the matrices of right or left eigenvectors of
A.
Arguments
=========
SIDE (input) CHARACTER*1
= 'R': compute right eigenvectors only;
= 'L': compute left eigenvectors only;
= 'B': compute both right and left eigenvectors.
HOWMNY (input) CHARACTER*1
= 'A': compute all right and/or left eigenvectors;
= 'B': compute all right and/or left eigenvectors, and
backtransform them using the input matrices supplied
in VR and/or VL;
= 'S': compute selected right and/or left eigenvectors,
specified by the logical array SELECT.
SELECT (input) LOGICAL array, dimension (N)
If HOWMNY='S', SELECT specifies the eigenvectors to be
computed.
If HOWMNY='A' or 'B', SELECT is not referenced.
To select the eigenvector corresponding to the j-th
eigenvalue, SELECT(j) must be set to .TRUE..
N (input) INTEGER
The order of the matrices A and B. N >= 0.
A (input) COMPLEX*16 array, dimension (LDA,N)
The upper triangular matrix A.
LDA (input) INTEGER
The leading dimension of array A. LDA >= max(1,N).
B (input) COMPLEX*16 array, dimension (LDB,N)
The upper triangular matrix B. B must have real diagonal
elements.
LDB (input) INTEGER
The leading dimension of array B. LDB >= max(1,N).
VL (input/output) 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 left Schur vectors returned by ZHGEQZ).
On exit, if SIDE = 'L' or 'B', VL contains:
if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B);
if HOWMNY = 'B', the matrix Q*Y;
if HOWMNY = 'S', the left eigenvectors of (A,B) specified by
SELECT, stored consecutively in the columns of
VL, in the same order as their eigenvalues.
If SIDE = 'R', VL is not referenced.
LDVL (input) INTEGER
The leading dimension of array VL.
LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
VR (input/output) 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 Z
of right Schur vectors returned by ZHGEQZ).
On exit, if SIDE = 'R' or 'B', VR contains:
if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B);
if HOWMNY = 'B', the matrix Z*X;
if HOWMNY = 'S', the right eigenvectors of (A,B) specified by
SELECT, stored consecutively in the columns of
VR, in the same order as their eigenvalues.
If SIDE = 'L', VR is not referenced.
LDVR (input) INTEGER
The leading dimension of the array VR.
LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
MM (input) INTEGER
The leading dimension of the array VR.
LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
MM (input) INTEGER
The number of columns in the arrays VL and/or VR. MM >= M.
M (output) 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.
WORK (workspace) COMPLEX*16 array, dimension (2*N)
RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
INFO (output) INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal value.
=====================================================================
Decode and Test the input parameters
Parameter adjustments
Function Body */
/* Table of constant values */
static doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__1 = 1;
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
vr_offset, i__1, i__2, i__3, i__4, i__5;
doublereal d__1, d__2, d__3, d__4, d__5, d__6;
doublecomplex z__1, z__2, z__3, z__4;
/* Builtin functions */
double d_imag(doublecomplex *);
void d_cnjg(doublecomplex *, doublecomplex *), z_div(doublecomplex *,
doublecomplex *, doublecomplex *);
/* Local variables */
static integer ibeg, ieig, iend;
static doublereal dmin__;
static integer isrc;
static doublereal temp;
static doublecomplex suma, sumb;
static doublereal xmax;
static doublecomplex d;
static integer i, j;
static doublereal scale;
static logical ilall;
static integer iside;
static doublereal sbeta;
extern logical lsame_(char *, char *);
static doublereal small;
static logical compl;
static doublereal anorm, bnorm;
static logical compr;
extern /* Subroutine */ int zgemv_(char *, integer *, integer *,
doublecomplex *, doublecomplex *, integer *, doublecomplex *,
integer *, doublecomplex *, doublecomplex *, integer *);
static doublecomplex ca, cb;
extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
static logical ilbbad;
static doublereal acoefa;
static integer je;
static doublereal bcoefa, acoeff;
static doublecomplex bcoeff;
static logical ilback;
static integer im;
static doublereal ascale, bscale;
extern doublereal dlamch_(char *);
static integer jr;
static doublecomplex salpha;
static doublereal safmin;
extern /* Subroutine */ int xerbla_(char *, integer *);
static doublereal bignum;
static logical ilcomp;
static integer ihwmny;
static doublereal big;
static logical lsa, lsb;
static doublereal ulp;
static doublecomplex sum;
#define SELECT(I) select[(I)-1]
#define WORK(I) work[(I)-1]
#define RWORK(I) rwork[(I)-1]
#define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
#define B(I,J) b[(I)-1 + ((J)-1)* ( *ldb)]
#define VL(I,J) vl[(I)-1 + ((J)-1)* ( *ldvl)]
#define VR(I,J) vr[(I)-1 + ((J)-1)* ( *ldvr)]
if (lsame_(howmny, "A")) {
ihwmny = 1;
ilall = TRUE_;
ilback = FALSE_;
} else if (lsame_(howmny, "S")) {
ihwmny = 2;
ilall = FALSE_;
ilback = FALSE_;
} else if (lsame_(howmny, "B") || lsame_(howmny, "T")) {
ihwmny = 3;
ilall = TRUE_;
ilback = TRUE_;
} else {
ihwmny = -1;
}
if (lsame_(side, "R")) {
iside = 1;
compl = FALSE_;
compr = TRUE_;
} else if (lsame_(side, "L")) {
iside = 2;
compl = TRUE_;
compr = FALSE_;
} else if (lsame_(side, "B")) {
iside = 3;
compl = TRUE_;
compr = TRUE_;
} else {
iside = -1;
}
/* Count the number of eigenvectors */
if (! ilall) {
im = 0;
i__1 = *n;
for (j = 1; j <= *n; ++j) {
if (SELECT(j)) {
++im;
}
/* L10: */
}
} else {
im = *n;
}
/* Check diagonal of B */
ilbbad = FALSE_;
i__1 = *n;
for (j = 1; j <= *n; ++j) {
if (d_imag(&B(j,j)) != 0.) {
ilbbad = TRUE_;
}
/* L20: */
}
*info = 0;
if (iside < 0) {
*info = -1;
} else if (ihwmny < 0) {
*info = -2;
} else if (*n < 0) {
*info = -4;
} else if (*lda < max(1,*n)) {
*info = -6;
} else if (ilbbad) {
*info = -7;
} else if (*ldb < max(1,*n)) {
*info = -8;
} else if (compl && *ldvl < *n || *ldvl < 1) {
*info = -10;
} else if (compr && *ldvr < *n || *ldvr < 1) {
*info = -12;
} else if (*mm < im) {
*info = -13;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("ZTGEVC", &i__1);
return 0;
}
/* Quick return if possible */
*m = im;
if (*n == 0) {
return 0;
}
/* Machine Constants */
safmin = dlamch_("Safe minimum");
big = 1. / safmin;
dlabad_(&safmin, &big);
ulp = dlamch_("Epsilon") * dlamch_("Base");
small = safmin * *n / ulp;
big = 1. / small;
bignum = 1. / (safmin * *n);
/* Compute the 1-norm of each column of the strictly upper triangular
part of A and B to check for possible overflow in the triangular
solver. */
i__1 = a_dim1 + 1;
anorm = (d__1 = A(1,1).r, abs(d__1)) + (d__2 = d_imag(&A(1,1)),
abs(d__2));
i__1 = b_dim1 + 1;
bnorm = (d__1 = B(1,1).r, abs(d__1)) + (d__2 = d_imag(&B(1,1)),
abs(d__2));
RWORK(1) = 0.;
RWORK(*n + 1) = 0.;
i__1 = *n;
for (j = 2; j <= *n; ++j) {
RWORK(j) = 0.;
RWORK(*n + j) = 0.;
i__2 = j - 1;
for (i = 1; i <= j-1; ++i) {
i__3 = i + j * a_dim1;
RWORK(j) += (d__1 = A(i,j).r, abs(d__1)) + (d__2 = d_imag(&A(i,j)), abs(d__2));
i__3 = i + j * b_dim1;
RWORK(*n + j) += (d__1 = B(i,j).r, abs(d__1)) + (d__2 = d_imag(&
B(i,j)), abs(d__2));
/* L30: */
}
/* Computing MAX */
i__2 = j + j * a_dim1;
d__3 = anorm, d__4 = RWORK(j) + ((d__1 = A(j,j).r, abs(d__1)) + (
d__2 = d_imag(&A(j,j)), abs(d__2)));
anorm = max(d__3,d__4);
/* Computing MAX */
i__2 = j + j * b_dim1;
d__3 = bnorm, d__4 = RWORK(*n + j) + ((d__1 = B(j,j).r, abs(d__1)) +
(d__2 = d_imag(&B(j,j)), abs(d__2)));
bnorm = max(d__3,d__4);
/* L40: */
}
ascale = 1. / max(anorm,safmin);
bscale = 1. / max(bnorm,safmin);
/* Left eigenvectors */
if (compl) {
ieig = 0;
/* Main loop over eigenvalues */
i__1 = *n;
for (je = 1; je <= *n; ++je) {
if (ilall) {
ilcomp = TRUE_;
} else {
ilcomp = SELECT(je);
}
if (ilcomp) {
++ieig;
i__2 = je + je * a_dim1;
i__3 = je + je * b_dim1;
if ((d__1 = A(je,je).r, abs(d__1)) + (d__2 = d_imag(&A(je,je)), abs(d__2)) <= safmin && (d__3 = B(je,je).r,
abs(d__3)) <= safmin) {
/* Singular matrix pencil -- return unit e
igenvector */
i__2 = *n;
for (jr = 1; jr <= *n; ++jr) {
i__3 = jr + ieig * vl_dim1;
VL(jr,ieig).r = 0., VL(jr,ieig).i = 0.;
/* L50: */
}
i__2 = ieig + ieig * vl_dim1;
VL(ieig,ieig).r = 1., VL(ieig,ieig).i = 0.;
goto L140;
}
/* Non-singular eigenvalue:
Compute coefficients a and b in
H
y ( a A - b B ) = 0
Computing MAX */
i__2 = je + je * a_dim1;
i__3 = je + je * b_dim1;
d__4 = ((d__1 = A(je,je).r, abs(d__1)) + (d__2 = d_imag(&A(je,je)), abs(d__2))) * ascale, d__5 = (d__3 =
B(je,je).r, abs(d__3)) * bscale, d__4 = max(d__4,d__5);
temp = 1. / max(d__4,safmin);
i__2 = je + je * a_dim1;
z__2.r = temp * A(je,je).r, z__2.i = temp * A(je,je).i;
z__1.r = ascale * z__2.r, z__1.i = ascale * z__2.i;
salpha.r = z__1.r, salpha.i = z__1.i;
i__2 = je + je * b_dim1;
sbeta = temp * B(je,je).r * bscale;
acoeff = sbeta * ascale;
z__1.r = bscale * salpha.r, z__1.i = bscale * salpha.i;
bcoeff.r = z__1.r, bcoeff.i = z__1.i;
/* Scale to avoid underflow */
lsa = abs(sbeta) >= safmin && abs(acoeff) < small;
lsb = (d__1 = salpha.r, abs(d__1)) + (d__2 = d_imag(&salpha),
abs(d__2)) >= safmin && (d__3 = bcoeff.r, abs(d__3))
+ (d__4 = d_imag(&bcoeff), abs(d__4)) < small;
scale = 1.;
if (lsa) {
scale = small / abs(sbeta) * min(anorm,big);
}
if (lsb) {
/* Computing MAX */
d__3 = scale, d__4 = small / ((d__1 = salpha.r, abs(d__1))
+ (d__2 = d_imag(&salpha), abs(d__2))) * min(
bnorm,big);
scale = max(d__3,d__4);
}
if (lsa || lsb) {
/* Computing MIN
Computing MAX */
d__5 = 1., d__6 = abs(acoeff), d__5 = max(d__5,d__6),
d__6 = (d__1 = bcoeff.r, abs(d__1)) + (d__2 =
d_imag(&bcoeff), abs(d__2));
d__3 = scale, d__4 = 1. / (safmin * max(d__5,d__6));
scale = min(d__3,d__4);
if (lsa) {
acoeff = ascale * (scale * sbeta);
} else {
acoeff = scale * acoeff;
}
if (lsb) {
z__2.r = scale * salpha.r, z__2.i = scale * salpha.i;
z__1.r = bscale * z__2.r, z__1.i = bscale * z__2.i;
bcoeff.r = z__1.r, bcoeff.i = z__1.i;
} else {
z__1.r = scale * bcoeff.r, z__1.i = scale * bcoeff.i;
bcoeff.r = z__1.r, bcoeff.i = z__1.i;
}
}
acoefa = abs(acoeff);
bcoefa = (d__1 = bcoeff.r, abs(d__1)) + (d__2 = d_imag(&
bcoeff), abs(d__2));
xmax = 1.;
i__2 = *n;
for (jr = 1; jr <= *n; ++jr) {
i__3 = jr;
WORK(jr).r = 0., WORK(jr).i = 0.;
/* L60: */
}
i__2 = je;
WORK(je).r = 1., WORK(je).i = 0.;
/* Computing MAX */
d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm,
d__1 = max(d__1,d__2);
dmin__ = max(d__1,safmin);
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?