ztgevc.c
来自「算断裂的」· C语言 代码 · 共 948 行 · 第 1/2 页
C
948 行
/* H
Triangular solve of (a A - b B) y = 0
H
(rowwise in (a A - b B) , or columnwise in a
A - b B) */
i__2 = *n;
for (j = je + 1; j <= *n; ++j) {
/* Compute
j-1
SUM = sum conjg( a*A(k,j) - b*B(k,j) )
*x(k)
k=je
(Scale if necessary) */
temp = 1. / xmax;
if (acoefa * RWORK(j) + bcoefa * RWORK(*n + j) > bignum *
temp) {
i__3 = j - 1;
for (jr = je; jr <= j-1; ++jr) {
i__4 = jr;
i__5 = jr;
z__1.r = temp * WORK(jr).r, z__1.i = temp *
WORK(jr).i;
WORK(jr).r = z__1.r, WORK(jr).i = z__1.i;
/* L70: */
}
xmax = 1.;
}
suma.r = 0., suma.i = 0.;
sumb.r = 0., sumb.i = 0.;
i__3 = j - 1;
for (jr = je; jr <= j-1; ++jr) {
d_cnjg(&z__3, &A(jr,j));
i__4 = jr;
z__2.r = z__3.r * WORK(jr).r - z__3.i * WORK(jr)
.i, z__2.i = z__3.r * WORK(jr).i + z__3.i *
WORK(jr).r;
z__1.r = suma.r + z__2.r, z__1.i = suma.i + z__2.i;
suma.r = z__1.r, suma.i = z__1.i;
d_cnjg(&z__3, &B(jr,j));
i__4 = jr;
z__2.r = z__3.r * WORK(jr).r - z__3.i * WORK(jr)
.i, z__2.i = z__3.r * WORK(jr).i + z__3.i *
WORK(jr).r;
z__1.r = sumb.r + z__2.r, z__1.i = sumb.i + z__2.i;
sumb.r = z__1.r, sumb.i = z__1.i;
/* L80: */
}
z__2.r = acoeff * suma.r, z__2.i = acoeff * suma.i;
d_cnjg(&z__4, &bcoeff);
z__3.r = z__4.r * sumb.r - z__4.i * sumb.i, z__3.i =
z__4.r * sumb.i + z__4.i * sumb.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
sum.r = z__1.r, sum.i = z__1.i;
/* Form x(j) = - SUM / conjg( a*A(j,j) - b
*B(j,j) )
with scaling and perturbation of the de
nominator */
i__3 = j + j * a_dim1;
z__3.r = acoeff * A(j,j).r, z__3.i = acoeff * A(j,j).i;
i__4 = j + j * b_dim1;
z__4.r = bcoeff.r * B(j,j).r - bcoeff.i * B(j,j).i,
z__4.i = bcoeff.r * B(j,j).i + bcoeff.i * B(j,j)
.r;
z__2.r = z__3.r - z__4.r, z__2.i = z__3.i - z__4.i;
d_cnjg(&z__1, &z__2);
d.r = z__1.r, d.i = z__1.i;
if ((d__1 = d.r, abs(d__1)) + (d__2 = d_imag(&d), abs(
d__2)) <= dmin__) {
z__1.r = dmin__, z__1.i = 0.;
d.r = z__1.r, d.i = z__1.i;
}
if ((d__1 = d.r, abs(d__1)) + (d__2 = d_imag(&d), abs(
d__2)) < 1.) {
if ((d__1 = sum.r, abs(d__1)) + (d__2 = d_imag(&sum),
abs(d__2)) >= bignum * ((d__3 = d.r, abs(d__3)
) + (d__4 = d_imag(&d), abs(d__4)))) {
temp = 1. / ((d__1 = sum.r, abs(d__1)) + (d__2 =
d_imag(&sum), abs(d__2)));
i__3 = j - 1;
for (jr = je; jr <= j-1; ++jr) {
i__4 = jr;
i__5 = jr;
z__1.r = temp * WORK(jr).r, z__1.i = temp *
WORK(jr).i;
WORK(jr).r = z__1.r, WORK(jr).i = z__1.i;
/* L90: */
}
xmax = temp * xmax;
z__1.r = temp * sum.r, z__1.i = temp * sum.i;
sum.r = z__1.r, sum.i = z__1.i;
}
}
i__3 = j;
z__2.r = -sum.r, z__2.i = -sum.i;
z_div(&z__1, &z__2, &d);
WORK(j).r = z__1.r, WORK(j).i = z__1.i;
/* Computing MAX */
i__3 = j;
d__3 = xmax, d__4 = (d__1 = WORK(j).r, abs(d__1)) + (
d__2 = d_imag(&WORK(j)), abs(d__2));
xmax = max(d__3,d__4);
/* L100: */
}
/* Back transform eigenvector if HOWMNY='B'. */
if (ilback) {
i__2 = *n + 1 - je;
zgemv_("N", n, &i__2, &c_b2, &VL(1,je), ldvl,
&WORK(je), &c__1, &c_b1, &WORK(*n + 1), &c__1)
;
isrc = 2;
ibeg = 1;
} else {
isrc = 1;
ibeg = je;
}
/* Copy and scale eigenvector into column of VL
*/
xmax = 0.;
i__2 = *n;
for (jr = ibeg; jr <= *n; ++jr) {
/* Computing MAX */
i__3 = (isrc - 1) * *n + jr;
d__3 = xmax, d__4 = (d__1 = WORK((isrc-1)**n+jr).r, abs(d__1)) + (
d__2 = d_imag(&WORK((isrc - 1) * *n + jr)), abs(
d__2));
xmax = max(d__3,d__4);
/* L110: */
}
if (xmax > safmin) {
temp = 1. / xmax;
i__2 = *n;
for (jr = ibeg; jr <= *n; ++jr) {
i__3 = jr + ieig * vl_dim1;
i__4 = (isrc - 1) * *n + jr;
z__1.r = temp * WORK((isrc-1)**n+jr).r, z__1.i = temp * WORK(
(isrc-1)**n+jr).i;
VL(jr,ieig).r = z__1.r, VL(jr,ieig).i = z__1.i;
/* L120: */
}
} else {
ibeg = *n + 1;
}
i__2 = ibeg - 1;
for (jr = 1; jr <= ibeg-1; ++jr) {
i__3 = jr + ieig * vl_dim1;
VL(jr,ieig).r = 0., VL(jr,ieig).i = 0.;
/* L130: */
}
}
L140:
;
}
}
/* Right eigenvectors */
if (compr) {
ieig = im + 1;
/* Main loop over eigenvalues */
for (je = *n; je >= 1; --je) {
if (ilall) {
ilcomp = TRUE_;
} else {
ilcomp = SELECT(je);
}
if (ilcomp) {
--ieig;
i__1 = je + je * a_dim1;
i__2 = 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__1 = *n;
for (jr = 1; jr <= *n; ++jr) {
i__2 = jr + ieig * vr_dim1;
VR(jr,ieig).r = 0., VR(jr,ieig).i = 0.;
/* L150: */
}
i__1 = ieig + ieig * vr_dim1;
VR(ieig,ieig).r = 1., VR(ieig,ieig).i = 0.;
goto L250;
}
/* Non-singular eigenvalue:
Compute coefficients a and b in
( a A - b B ) x = 0
Computing MAX */
i__1 = je + je * a_dim1;
i__2 = 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__1 = 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__1 = 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__1 = *n;
for (jr = 1; jr <= *n; ++jr) {
i__2 = jr;
WORK(jr).r = 0., WORK(jr).i = 0.;
/* L160: */
}
i__1 = 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);
/* Triangular solve of (a A - b B) x = 0 (colum
nwise)
WORK(1:j-1) contains sums w,
WORK(j+1:JE) contains x */
i__1 = je - 1;
for (jr = 1; jr <= je-1; ++jr) {
i__2 = jr;
i__3 = jr + je * a_dim1;
z__2.r = acoeff * A(jr,je).r, z__2.i = acoeff * A(jr,je).i;
i__4 = jr + je * b_dim1;
z__3.r = bcoeff.r * B(jr,je).r - bcoeff.i * B(jr,je).i,
z__3.i = bcoeff.r * B(jr,je).i + bcoeff.i * B(jr,je)
.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
WORK(jr).r = z__1.r, WORK(jr).i = z__1.i;
/* L170: */
}
i__1 = je;
WORK(je).r = 1., WORK(je).i = 0.;
for (j = je - 1; j >= 1; --j) {
/* Form x(j) := - w(j) / d
with scaling and perturbation of the de
nominator */
i__1 = j + j * a_dim1;
z__2.r = acoeff * A(j,j).r, z__2.i = acoeff * A(j,j).i;
i__2 = j + j * b_dim1;
z__3.r = bcoeff.r * B(j,j).r - bcoeff.i * B(j,j).i,
z__3.i = bcoeff.r * B(j,j).i + bcoeff.i * B(j,j)
.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
d.r = z__1.r, d.i = z__1.i;
if ((d__1 = d.r, abs(d__1)) + (d__2 = d_imag(&d), abs(
d__2)) <= dmin__) {
z__1.r = dmin__, z__1.i = 0.;
d.r = z__1.r, d.i = z__1.i;
}
if ((d__1 = d.r, abs(d__1)) + (d__2 = d_imag(&d), abs(
d__2)) < 1.) {
i__1 = j;
if ((d__1 = WORK(j).r, abs(d__1)) + (d__2 = d_imag(
&WORK(j)), abs(d__2)) >= bignum * ((d__3 =
d.r, abs(d__3)) + (d__4 = d_imag(&d), abs(
d__4)))) {
i__1 = j;
temp = 1. / ((d__1 = WORK(j).r, abs(d__1)) + (
d__2 = d_imag(&WORK(j)), abs(d__2)));
i__1 = je;
for (jr = 1; jr <= je; ++jr) {
i__2 = jr;
i__3 = jr;
z__1.r = temp * WORK(jr).r, z__1.i = temp *
WORK(jr).i;
WORK(jr).r = z__1.r, WORK(jr).i = z__1.i;
/* L180: */
}
}
}
i__1 = j;
i__2 = j;
z__2.r = -WORK(j).r, z__2.i = -WORK(j).i;
z_div(&z__1, &z__2, &d);
WORK(j).r = z__1.r, WORK(j).i = z__1.i;
if (j > 1) {
/* w = w + x(j)*(a A(*,j) - b B(*,j
) ) with scaling */
i__1 = j;
if ((d__1 = WORK(j).r, abs(d__1)) + (d__2 = d_imag(
&WORK(j)), abs(d__2)) > 1.) {
i__1 = j;
temp = 1. / ((d__1 = WORK(j).r, abs(d__1)) + (
d__2 = d_imag(&WORK(j)), abs(d__2)));
if (acoefa * RWORK(j) + bcoefa * RWORK(*n + j) >=
bignum * temp) {
i__1 = je;
for (jr = 1; jr <= je; ++jr) {
i__2 = jr;
i__3 = jr;
z__1.r = temp * WORK(jr).r, z__1.i =
temp * WORK(jr).i;
WORK(jr).r = z__1.r, WORK(jr).i =
z__1.i;
/* L190: */
}
}
}
i__1 = j;
z__1.r = acoeff * WORK(j).r, z__1.i = acoeff *
WORK(j).i;
ca.r = z__1.r, ca.i = z__1.i;
i__1 = j;
z__1.r = bcoeff.r * WORK(j).r - bcoeff.i * WORK(
j).i, z__1.i = bcoeff.r * WORK(j).i +
bcoeff.i * WORK(j).r;
cb.r = z__1.r, cb.i = z__1.i;
i__1 = j - 1;
for (jr = 1; jr <= j-1; ++jr) {
i__2 = jr;
i__3 = jr;
i__4 = jr + j * a_dim1;
z__3.r = ca.r * A(jr,j).r - ca.i * A(jr,j).i,
z__3.i = ca.r * A(jr,j).i + ca.i * A(jr,j)
.r;
z__2.r = WORK(jr).r + z__3.r, z__2.i = WORK(
jr).i + z__3.i;
i__5 = jr + j * b_dim1;
z__4.r = cb.r * B(jr,j).r - cb.i * B(jr,j).i,
z__4.i = cb.r * B(jr,j).i + cb.i * B(jr,j)
.r;
z__1.r = z__2.r - z__4.r, z__1.i = z__2.i -
z__4.i;
WORK(jr).r = z__1.r, WORK(jr).i = z__1.i;
/* L200: */
}
}
/* L210: */
}
/* Back transform eigenvector if HOWMNY='B'. */
if (ilback) {
zgemv_("N", n, &je, &c_b2, &VR(1,1), ldvr, &WORK(1),
&c__1, &c_b1, &WORK(*n + 1), &c__1);
isrc = 2;
iend = *n;
} else {
isrc = 1;
iend = je;
}
/* Copy and scale eigenvector into column of VR
*/
xmax = 0.;
i__1 = iend;
for (jr = 1; jr <= iend; ++jr) {
/* Computing MAX */
i__2 = (isrc - 1) * *n + jr;
d__3 = xmax, d__4 = (d__1 = WORK((isrc-1)**n+jr).r, abs(d__1)) + (
d__2 = d_imag(&WORK((isrc - 1) * *n + jr)), abs(
d__2));
xmax = max(d__3,d__4);
/* L220: */
}
if (xmax > safmin) {
temp = 1. / xmax;
i__1 = iend;
for (jr = 1; jr <= iend; ++jr) {
i__2 = jr + ieig * vr_dim1;
i__3 = (isrc - 1) * *n + jr;
z__1.r = temp * WORK((isrc-1)**n+jr).r, z__1.i = temp * WORK(
(isrc-1)**n+jr).i;
VR(jr,ieig).r = z__1.r, VR(jr,ieig).i = z__1.i;
/* L230: */
}
} else {
iend = 0;
}
i__1 = *n;
for (jr = iend + 1; jr <= *n; ++jr) {
i__2 = jr + ieig * vr_dim1;
VR(jr,ieig).r = 0., VR(jr,ieig).i = 0.;
/* L240: */
}
}
L250:
;
}
}
return 0;
/* End of ZTGEVC */
} /* ztgevc_ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?