📄 dhgeqz.c
字号:
tempr = sqr * szr - sqi * szi;
tempi = sqr * szi + sqi * szr;
b1r = cq * cz * b11 + tempr * b22;
b1i = tempi * b22;
b1a = dlapy2_(&b1r, &b1i);
b2r = cq * cz * b22 + tempr * b11;
b2i = -tempi * b11;
b2a = dlapy2_(&b2r, &b2i);
/* Normalize so beta > 0, and Im( alpha1 ) > 0 */
beta[ilast - 1] = b1a;
beta[ilast] = b2a;
alphar[ilast - 1] = wr * b1a * s1inv;
alphai[ilast - 1] = wi * b1a * s1inv;
alphar[ilast] = wr * b2a * s1inv;
alphai[ilast] = -(wi * b2a) * s1inv;
/* Step 3: Go to next block -- exit if finished. */
ilast = ifirst - 1;
if (ilast < *ilo) {
goto L380;
}
/* Reset counters */
iiter = 0;
eshift = 0.;
if (! ilschr) {
ilastm = ilast;
if (ifrstm > ilast) {
ifrstm = *ilo;
}
}
} else {
/* Usual case: 3x3 or larger block, using Francis implicit */
/* double-shift */
/* 2 */
/* Eigenvalue equation is w - c w + d = 0, */
/* -1 2 -1 */
/* so compute 1st column of (A B ) - c A B + d */
/* using the formula in QZIT (from EISPACK) */
/* We assume that the block is at least 3x3 */
ad11 = ascale * a[ilast - 1 + (ilast - 1) * a_dim1] / (bscale * b[ilast - 1 + (ilast - 1) * b_dim1]);
ad21 = ascale * a[ilast + (ilast - 1) * a_dim1] / (bscale * b[ilast - 1 + (ilast - 1) * b_dim1]);
ad12 = ascale * a[ilast - 1 + ilast * a_dim1] / (bscale * b[ilast + ilast * b_dim1]);
ad22 = ascale * a[ilast + ilast * a_dim1] / (bscale * b[ilast + ilast * b_dim1]);
u12 = b[ilast - 1 + ilast * b_dim1] / b[ilast + ilast * b_dim1];
ad11l = ascale * a[ifirst + ifirst * a_dim1] / (bscale * b[ifirst + ifirst * b_dim1]);
ad21l = ascale * a[ifirst + 1 + ifirst * a_dim1] / (bscale * b[ifirst + ifirst * b_dim1]);
ad12l = ascale * a[ifirst + (ifirst + 1) * a_dim1] / (bscale * b[ifirst + 1 + (ifirst + 1) * b_dim1]);
ad22l = ascale * a[ifirst + 1 + (ifirst + 1) * a_dim1] / (bscale * b[ifirst + 1 + (ifirst + 1) * b_dim1]);
ad32l = ascale * a[ifirst + 2 + (ifirst + 1) * a_dim1] / (bscale * b[ifirst + 1 + (ifirst + 1) * b_dim1]);
u12l = b[ifirst + (ifirst + 1) * b_dim1] / b[ifirst + 1 + (ifirst + 1) * b_dim1];
v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12 * ad11l + (ad12l - ad11l * u12l) * ad21l;
v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 - ad11l) + ad21 * u12) * ad21l;
v[2] = ad32l * ad21l;
istart = ifirst;
dlarfg_(&c__3, v, &v[1], &c__1, &tau);
v[0] = 1.;
/* Sweep */
for (j = istart; j < ilast - 1; ++j) {
/* All but last elements: use 3x3 Householder transforms. */
/* Zero (j-1)st column of A */
if (j > istart) {
v[0] = a[j + (j - 1) * a_dim1];
v[1] = a[j + 1 + (j - 1) * a_dim1];
v[2] = a[j + 2 + (j - 1) * a_dim1];
dlarfg_(&c__3, &a[j + (j - 1) * a_dim1], &v[1], &c__1, &tau);
v[0] = 1.;
a[j + 1 + (j - 1) * a_dim1] = 0.;
a[j + 2 + (j - 1) * a_dim1] = 0.;
}
for (jc = j; jc <= ilastm; ++jc) {
temp = tau * (a[j + jc * a_dim1] + v[1] * a[j + 1 + jc * a_dim1] + v[2] * a[j + 2 + jc * a_dim1]);
a[j + jc * a_dim1] -= temp;
a[j + 1 + jc * a_dim1] -= temp * v[1];
a[j + 2 + jc * a_dim1] -= temp * v[2];
temp2 = tau * (b[j + jc * b_dim1] + v[1] * b[j + 1 + jc * b_dim1] + v[2] * b[j + 2 + jc * b_dim1]);
b[j + jc * b_dim1] -= temp2;
b[j + 1 + jc * b_dim1] -= temp2 * v[1];
b[j + 2 + jc * b_dim1] -= temp2 * v[2];
}
if (ilq) {
for (jr = 1; jr <= *n; ++jr) {
temp = tau * (q[jr + j * q_dim1] + v[1] * q[jr + (j + 1) * q_dim1] + v[2] * q[jr + (j + 2) * q_dim1]);
q[jr + j * q_dim1] -= temp;
q[jr + (j + 1) * q_dim1] -= temp * v[1];
q[jr + (j + 2) * q_dim1] -= temp * v[2];
}
}
/* Zero j-th column of B (see DLAGBC for details) */
/* Swap rows to pivot */
ilpivt = FALSE_;
temp = max(abs(b[j+1 + (j+1) * b_dim1]), abs(b[j+1 + (j+2) * b_dim1]));
temp2 = max(abs(b[j+2 + (j+1) * b_dim1]), abs(b[j+2 + (j+2) * b_dim1]));
if (max(temp,temp2) < safmin) {
scale = 0.;
u1 = 1.;
u2 = 0.;
goto L250;
} else if (temp >= temp2) {
w11 = b[j + 1 + (j + 1) * b_dim1];
w21 = b[j + 2 + (j + 1) * b_dim1];
w12 = b[j + 1 + (j + 2) * b_dim1];
w22 = b[j + 2 + (j + 2) * b_dim1];
u1 = b[j + 1 + j * b_dim1];
u2 = b[j + 2 + j * b_dim1];
} else {
w21 = b[j + 1 + (j + 1) * b_dim1];
w11 = b[j + 2 + (j + 1) * b_dim1];
w22 = b[j + 1 + (j + 2) * b_dim1];
w12 = b[j + 2 + (j + 2) * b_dim1];
u2 = b[j + 1 + j * b_dim1];
u1 = b[j + 2 + j * b_dim1];
}
/* Swap columns if nec. */
if (abs(w12) > abs(w11)) {
ilpivt = TRUE_;
temp = w12;
temp2 = w22;
w12 = w11;
w22 = w21;
w11 = temp;
w21 = temp2;
}
/* LU-factor */
temp = w21 / w11;
u2 -= temp * u1;
w22 -= temp * w12;
w21 = 0.;
/* Compute SCALE */
scale = 1.;
if (abs(w22) < safmin) {
scale = 0.;
u2 = 1.;
u1 = -w12 / w11;
goto L250;
}
if (abs(w22) < abs(u2)) {
scale = abs(w22 / u2);
}
if (abs(w11) < abs(u1)) {
scale = min(scale, abs(w11/u1));
}
/* Solve */
u2 = scale * u2 / w22;
u1 = (scale * u1 - w12 * u2) / w11;
L250:
if (ilpivt) {
temp = u2;
u2 = u1;
u1 = temp;
}
/* Compute Householder Vector */
t = sqrt(scale * scale + u1 * u1 + u2 * u2);
tau = scale / t + 1.;
vs = -1. / (scale + t);
v[0] = 1.;
v[1] = vs * u1;
v[2] = vs * u2;
/* Apply transformations from the right. */
for (jr = ifrstm; jr <= min(j+3, ilast); ++jr) {
temp = tau * (a[jr + j * a_dim1] + v[1] * a[jr + (j + 1) * a_dim1] + v[2] * a[jr + (j + 2) * a_dim1]);
a[jr + j * a_dim1] -= temp;
a[jr + (j + 1) * a_dim1] -= temp * v[1];
a[jr + (j + 2) * a_dim1] -= temp * v[2];
}
for (jr = ifrstm; jr <= j+2; ++jr) {
temp = tau * (b[jr + j * b_dim1] + v[1] * b[jr + (j + 1) * b_dim1] + v[2] * b[jr + (j + 2) * b_dim1]);
b[jr + j * b_dim1] -= temp;
b[jr + (j + 1) * b_dim1] -= temp * v[1];
b[jr + (j + 2) * b_dim1] -= temp * v[2];
}
if (ilz) {
for (jr = 1; jr <= *n; ++jr) {
temp = tau * (z[jr + j * z_dim1] + v[1] * z[jr + (j + 1) * z_dim1] + v[2] * z[jr + (j + 2) * z_dim1]);
z[jr + j * z_dim1] -= temp;
z[jr + (j + 1) * z_dim1] -= temp * v[1];
z[jr + (j + 2) * z_dim1] -= temp * v[2];
}
}
b[j + 1 + j * b_dim1] = 0.;
b[j + 2 + j * b_dim1] = 0.;
}
/* Last elements: Use Givens rotations */
/* Rotations from the left */
j = ilast - 1;
temp = a[j + (j - 1) * a_dim1];
dlartg_(&temp, &a[j + 1 + (j - 1) * a_dim1], &c, &s, &a[j + (j - 1) * a_dim1]);
a[j + 1 + (j - 1) * a_dim1] = 0.;
for (jc = j; jc <= ilastm; ++jc) {
temp = c * a[j + jc * a_dim1] + s * a[j + 1 + jc * a_dim1];
a[j + 1 + jc * a_dim1] = -s * a[j + jc * a_dim1] + c * a[j + 1 + jc * a_dim1];
a[j + jc * a_dim1] = temp;
temp2 = c * b[j + jc * b_dim1] + s * b[j + 1 + jc * b_dim1];
b[j + 1 + jc * b_dim1] = -s * b[j + jc * b_dim1] + c * b[j + 1 + jc * b_dim1];
b[j + jc * b_dim1] = temp2;
}
if (ilq) {
for (jr = 1; jr <= *n; ++jr) {
temp = c * q[jr + j * q_dim1] + s * q[jr + (j + 1) * q_dim1];
q[jr + (j + 1) * q_dim1] = -s * q[jr + j * q_dim1] + c * q[jr + (j + 1) * q_dim1];
q[jr + j * q_dim1] = temp;
}
}
/* Rotations from the right. */
temp = b[j + 1 + (j + 1) * b_dim1];
dlartg_(&temp, &b[j + 1 + j * b_dim1], &c, &s, &b[j + 1 + (j + 1) * b_dim1]);
b[j + 1 + j * b_dim1] = 0.;
for (jr = ifrstm; jr <= ilast; ++jr) {
temp = c * a[jr + (j + 1) * a_dim1] + s * a[jr + j * a_dim1];
a[jr + j * a_dim1] = -s * a[jr + (j + 1) * a_dim1] + c * a[jr + j * a_dim1];
a[jr + (j + 1) * a_dim1] = temp;
}
for (jr = ifrstm; jr < ilast; ++jr) {
temp = c * b[jr + (j + 1) * b_dim1] + s * b[jr + j * b_dim1];
b[jr + j * b_dim1] = -s * b[jr + (j + 1) * b_dim1] + c * b[jr + j * b_dim1];
b[jr + (j + 1) * b_dim1] = temp;
}
if (ilz) {
for (jr = 1; jr <= *n; ++jr) {
temp = c * z[jr + (j + 1) * z_dim1] + s * z[jr + j * z_dim1];
z[jr + j * z_dim1] = -s * z[jr + (j + 1) * z_dim1] + c * z[jr + j * z_dim1];
z[jr + (j + 1) * z_dim1] = temp;
}
}
/* End of Double-Shift code */
}
/* End of iteration loop */
L350:
;
}
/* Drop-through = non-convergence */
*info = ilast;
goto L420;
/* Successful completion of all QZ steps */
L380:
/* Set Eigenvalues 1:ILO-1 */
for (j = 1; j < *ilo; ++j) {
if (b[j + j * b_dim1] < 0.) {
if (ilschr) {
for (jr = 1; jr <= j; ++jr) {
a[jr + j * a_dim1] = -a[jr + j * a_dim1];
b[jr + j * b_dim1] = -b[jr + j * b_dim1];
}
} else {
a[j + j * a_dim1] = -a[j + j * a_dim1];
b[j + j * b_dim1] = -b[j + j * b_dim1];
}
if (ilz) {
for (jr = 1; jr <= *n; ++jr) {
z[jr + j * z_dim1] = -z[jr + j * z_dim1];
}
}
}
alphar[j] = a[j + j * a_dim1];
alphai[j] = 0.;
beta[j] = b[j + j * b_dim1];
}
/* Normal Termination */
*info = 0;
/* Exit (other than argument error) -- return optimal workspace size */
L420:
work[1] = (doublereal) (*n);
} /* dhgeqz_ */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -