📄 chgeqz.c
字号:
i__3 = b_subscr(ilast - 1, ilast - 1);
q__3.r = bscale * b[i__3].r, q__3.i = bscale * b[i__3].i;
c_div(&q__1, &q__2, &q__3);
ad21.r = q__1.r, ad21.i = q__1.i;
i__2 = a_subscr(ilast - 1, ilast);
q__2.r = ascale * a[i__2].r, q__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast, ilast);
q__3.r = bscale * b[i__3].r, q__3.i = bscale * b[i__3].i;
c_div(&q__1, &q__2, &q__3);
ad12.r = q__1.r, ad12.i = q__1.i;
i__2 = a_subscr(ilast, ilast);
q__2.r = ascale * a[i__2].r, q__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast, ilast);
q__3.r = bscale * b[i__3].r, q__3.i = bscale * b[i__3].i;
c_div(&q__1, &q__2, &q__3);
ad22.r = q__1.r, ad22.i = q__1.i;
q__2.r = u12.r * ad21.r - u12.i * ad21.i, q__2.i = u12.r * ad21.i
+ u12.i * ad21.r;
q__1.r = ad22.r - q__2.r, q__1.i = ad22.i - q__2.i;
abi22.r = q__1.r, abi22.i = q__1.i;
q__2.r = ad11.r + abi22.r, q__2.i = ad11.i + abi22.i;
q__1.r = q__2.r * .5f, q__1.i = q__2.i * .5f;
t.r = q__1.r, t.i = q__1.i;
pow_ci(&q__4, &t, &c__2);
q__5.r = ad12.r * ad21.r - ad12.i * ad21.i, q__5.i = ad12.r *
ad21.i + ad12.i * ad21.r;
q__3.r = q__4.r + q__5.r, q__3.i = q__4.i + q__5.i;
q__6.r = ad11.r * ad22.r - ad11.i * ad22.i, q__6.i = ad11.r *
ad22.i + ad11.i * ad22.r;
q__2.r = q__3.r - q__6.r, q__2.i = q__3.i - q__6.i;
c_sqrt(&q__1, &q__2);
rtdisc.r = q__1.r, rtdisc.i = q__1.i;
q__1.r = t.r - abi22.r, q__1.i = t.i - abi22.i;
q__2.r = t.r - abi22.r, q__2.i = t.i - abi22.i;
temp = q__1.r * rtdisc.r + r_imag(&q__2) * r_imag(&rtdisc);
if (temp <= 0.f) {
q__1.r = t.r + rtdisc.r, q__1.i = t.i + rtdisc.i;
shift.r = q__1.r, shift.i = q__1.i;
} else {
q__1.r = t.r - rtdisc.r, q__1.i = t.i - rtdisc.i;
shift.r = q__1.r, shift.i = q__1.i;
}
/* ------------------- Begin Timing Code ---------------------- */
opst += 116.f;
/* -------------------- End Timing Code ----------------------- */
} else {
/* Exceptional shift. Chosen for no particularly good reason. */
i__2 = a_subscr(ilast - 1, ilast);
q__4.r = ascale * a[i__2].r, q__4.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast - 1, ilast - 1);
q__5.r = bscale * b[i__3].r, q__5.i = bscale * b[i__3].i;
c_div(&q__3, &q__4, &q__5);
r_cnjg(&q__2, &q__3);
q__1.r = eshift.r + q__2.r, q__1.i = eshift.i + q__2.i;
eshift.r = q__1.r, eshift.i = q__1.i;
shift.r = eshift.r, shift.i = eshift.i;
/* ------------------- Begin Timing Code ---------------------- */
opst += 15.f;
/* -------------------- End Timing Code ----------------------- */
}
/* Now check for two consecutive small subdiagonals. */
i__2 = ifirst + 1;
for (j = ilast - 1; j >= i__2; --j) {
istart = j;
i__3 = a_subscr(j, j);
q__2.r = ascale * a[i__3].r, q__2.i = ascale * a[i__3].i;
i__4 = b_subscr(j, j);
q__4.r = bscale * b[i__4].r, q__4.i = bscale * b[i__4].i;
q__3.r = shift.r * q__4.r - shift.i * q__4.i, q__3.i = shift.r *
q__4.i + shift.i * q__4.r;
q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
temp = (r__1 = ctemp.r, dabs(r__1)) + (r__2 = r_imag(&ctemp),
dabs(r__2));
i__3 = a_subscr(j + 1, j);
temp2 = ascale * ((r__1 = a[i__3].r, dabs(r__1)) + (r__2 = r_imag(
&a_ref(j + 1, j)), dabs(r__2)));
tempr = dmax(temp,temp2);
if (tempr < 1.f && tempr != 0.f) {
temp /= tempr;
temp2 /= tempr;
}
i__3 = a_subscr(j, j - 1);
if (((r__1 = a[i__3].r, dabs(r__1)) + (r__2 = r_imag(&a_ref(j, j
- 1)), dabs(r__2))) * temp2 <= temp * atol) {
goto L90;
}
/* L80: */
}
istart = ifirst;
i__2 = a_subscr(ifirst, ifirst);
q__2.r = ascale * a[i__2].r, q__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ifirst, ifirst);
q__4.r = bscale * b[i__3].r, q__4.i = bscale * b[i__3].i;
q__3.r = shift.r * q__4.r - shift.i * q__4.i, q__3.i = shift.r *
q__4.i + shift.i * q__4.r;
q__1.r = q__2.r - q__3.r, q__1.i = q__2.i - q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
/* --------------------- Begin Timing Code ----------------------- */
opst += -6.f;
/* ---------------------- End Timing Code ------------------------ */
L90:
/* Do an implicit-shift QZ sweep.
Initial Q */
i__2 = a_subscr(istart + 1, istart);
q__1.r = ascale * a[i__2].r, q__1.i = ascale * a[i__2].i;
ctemp2.r = q__1.r, ctemp2.i = q__1.i;
/* --------------------- Begin Timing Code ----------------------- */
opst += (real) ((ilast - istart) * 18 + 2);
/* ---------------------- End Timing Code ------------------------ */
clartg_(&ctemp, &ctemp2, &c__, &s, &ctemp3);
/* Sweep */
i__2 = ilast - 1;
for (j = istart; j <= i__2; ++j) {
if (j > istart) {
i__3 = a_subscr(j, j - 1);
ctemp.r = a[i__3].r, ctemp.i = a[i__3].i;
clartg_(&ctemp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j -
1));
i__3 = a_subscr(j + 1, j - 1);
a[i__3].r = 0.f, a[i__3].i = 0.f;
}
i__3 = ilastm;
for (jc = j; jc <= i__3; ++jc) {
i__4 = a_subscr(j, jc);
q__2.r = c__ * a[i__4].r, q__2.i = c__ * a[i__4].i;
i__5 = a_subscr(j + 1, jc);
q__3.r = s.r * a[i__5].r - s.i * a[i__5].i, q__3.i = s.r * a[
i__5].i + s.i * a[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
i__4 = a_subscr(j + 1, jc);
r_cnjg(&q__4, &s);
q__3.r = -q__4.r, q__3.i = -q__4.i;
i__5 = a_subscr(j, jc);
q__2.r = q__3.r * a[i__5].r - q__3.i * a[i__5].i, q__2.i =
q__3.r * a[i__5].i + q__3.i * a[i__5].r;
i__6 = a_subscr(j + 1, jc);
q__5.r = c__ * a[i__6].r, q__5.i = c__ * a[i__6].i;
q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
a[i__4].r = q__1.r, a[i__4].i = q__1.i;
i__4 = a_subscr(j, jc);
a[i__4].r = ctemp.r, a[i__4].i = ctemp.i;
i__4 = b_subscr(j, jc);
q__2.r = c__ * b[i__4].r, q__2.i = c__ * b[i__4].i;
i__5 = b_subscr(j + 1, jc);
q__3.r = s.r * b[i__5].r - s.i * b[i__5].i, q__3.i = s.r * b[
i__5].i + s.i * b[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp2.r = q__1.r, ctemp2.i = q__1.i;
i__4 = b_subscr(j + 1, jc);
r_cnjg(&q__4, &s);
q__3.r = -q__4.r, q__3.i = -q__4.i;
i__5 = b_subscr(j, jc);
q__2.r = q__3.r * b[i__5].r - q__3.i * b[i__5].i, q__2.i =
q__3.r * b[i__5].i + q__3.i * b[i__5].r;
i__6 = b_subscr(j + 1, jc);
q__5.r = c__ * b[i__6].r, q__5.i = c__ * b[i__6].i;
q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
b[i__4].r = q__1.r, b[i__4].i = q__1.i;
i__4 = b_subscr(j, jc);
b[i__4].r = ctemp2.r, b[i__4].i = ctemp2.i;
/* L100: */
}
if (ilq) {
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
i__4 = q_subscr(jr, j);
q__2.r = c__ * q[i__4].r, q__2.i = c__ * q[i__4].i;
r_cnjg(&q__4, &s);
i__5 = q_subscr(jr, j + 1);
q__3.r = q__4.r * q[i__5].r - q__4.i * q[i__5].i, q__3.i =
q__4.r * q[i__5].i + q__4.i * q[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
i__4 = q_subscr(jr, j + 1);
q__3.r = -s.r, q__3.i = -s.i;
i__5 = q_subscr(jr, j);
q__2.r = q__3.r * q[i__5].r - q__3.i * q[i__5].i, q__2.i =
q__3.r * q[i__5].i + q__3.i * q[i__5].r;
i__6 = q_subscr(jr, j + 1);
q__4.r = c__ * q[i__6].r, q__4.i = c__ * q[i__6].i;
q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
q[i__4].r = q__1.r, q[i__4].i = q__1.i;
i__4 = q_subscr(jr, j);
q[i__4].r = ctemp.r, q[i__4].i = ctemp.i;
/* L110: */
}
}
i__3 = b_subscr(j + 1, j + 1);
ctemp.r = b[i__3].r, ctemp.i = b[i__3].i;
clartg_(&ctemp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
i__3 = b_subscr(j + 1, j);
b[i__3].r = 0.f, b[i__3].i = 0.f;
/* Computing MIN */
i__4 = j + 2;
i__3 = min(i__4,ilast);
for (jr = ifrstm; jr <= i__3; ++jr) {
i__4 = a_subscr(jr, j + 1);
q__2.r = c__ * a[i__4].r, q__2.i = c__ * a[i__4].i;
i__5 = a_subscr(jr, j);
q__3.r = s.r * a[i__5].r - s.i * a[i__5].i, q__3.i = s.r * a[
i__5].i + s.i * a[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
i__4 = a_subscr(jr, j);
r_cnjg(&q__4, &s);
q__3.r = -q__4.r, q__3.i = -q__4.i;
i__5 = a_subscr(jr, j + 1);
q__2.r = q__3.r * a[i__5].r - q__3.i * a[i__5].i, q__2.i =
q__3.r * a[i__5].i + q__3.i * a[i__5].r;
i__6 = a_subscr(jr, j);
q__5.r = c__ * a[i__6].r, q__5.i = c__ * a[i__6].i;
q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
a[i__4].r = q__1.r, a[i__4].i = q__1.i;
i__4 = a_subscr(jr, j + 1);
a[i__4].r = ctemp.r, a[i__4].i = ctemp.i;
/* L120: */
}
i__3 = j;
for (jr = ifrstm; jr <= i__3; ++jr) {
i__4 = b_subscr(jr, j + 1);
q__2.r = c__ * b[i__4].r, q__2.i = c__ * b[i__4].i;
i__5 = b_subscr(jr, j);
q__3.r = s.r * b[i__5].r - s.i * b[i__5].i, q__3.i = s.r * b[
i__5].i + s.i * b[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
i__4 = b_subscr(jr, j);
r_cnjg(&q__4, &s);
q__3.r = -q__4.r, q__3.i = -q__4.i;
i__5 = b_subscr(jr, j + 1);
q__2.r = q__3.r * b[i__5].r - q__3.i * b[i__5].i, q__2.i =
q__3.r * b[i__5].i + q__3.i * b[i__5].r;
i__6 = b_subscr(jr, j);
q__5.r = c__ * b[i__6].r, q__5.i = c__ * b[i__6].i;
q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
b[i__4].r = q__1.r, b[i__4].i = q__1.i;
i__4 = b_subscr(jr, j + 1);
b[i__4].r = ctemp.r, b[i__4].i = ctemp.i;
/* L130: */
}
if (ilz) {
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
i__4 = z___subscr(jr, j + 1);
q__2.r = c__ * z__[i__4].r, q__2.i = c__ * z__[i__4].i;
i__5 = z___subscr(jr, j);
q__3.r = s.r * z__[i__5].r - s.i * z__[i__5].i, q__3.i =
s.r * z__[i__5].i + s.i * z__[i__5].r;
q__1.r = q__2.r + q__3.r, q__1.i = q__2.i + q__3.i;
ctemp.r = q__1.r, ctemp.i = q__1.i;
i__4 = z___subscr(jr, j);
r_cnjg(&q__4, &s);
q__3.r = -q__4.r, q__3.i = -q__4.i;
i__5 = z___subscr(jr, j + 1);
q__2.r = q__3.r * z__[i__5].r - q__3.i * z__[i__5].i,
q__2.i = q__3.r * z__[i__5].i + q__3.i * z__[i__5]
.r;
i__6 = z___subscr(jr, j);
q__5.r = c__ * z__[i__6].r, q__5.i = c__ * z__[i__6].i;
q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
z__[i__4].r = q__1.r, z__[i__4].i = q__1.i;
i__4 = z___subscr(jr, j + 1);
z__[i__4].r = ctemp.r, z__[i__4].i = ctemp.i;
/* L140: */
}
}
/* L150: */
}
/* --------------------- Begin Timing Code ----------------------- */
opst += (real) (ilast - istart) * (real) ((ilastm - ifrstm) * 40 +
184 + (nq + nz) * 20) - 20;
/* ---------------------- End Timing Code ------------------------ */
L160:
/* --------------------- Begin Timing Code -----------------------
End of iteration -- add in "small" contributions. */
latime_1.ops += opst;
opst = 0.f;
/* ---------------------- End Timing Code ------------------------
L170: */
}
/* Drop-through = non-convergence */
L180:
*info = ilast;
/* ---------------------- Begin Timing Code ------------------------- */
latime_1.ops += opst;
opst = 0.f;
/* ----------------------- End Timing Code -------------------------- */
goto L210;
/* Successful completion of all QZ steps */
L190:
/* Set Eigenvalues 1:ILO-1 */
i__1 = *ilo - 1;
for (j = 1; j <= i__1; ++j) {
absb = c_abs(&b_ref(j, j));
if (absb > safmin) {
i__2 = b_subscr(j, j);
q__2.r = b[i__2].r / absb, q__2.i = b[i__2].i / absb;
r_cnjg(&q__1, &q__2);
signbc.r = q__1.r, signbc.i = q__1.i;
i__2 = b_subscr(j, j);
b[i__2].r = absb, b[i__2].i = 0.f;
if (ilschr) {
i__2 = j - 1;
cscal_(&i__2, &signbc, &b_ref(1, j), &c__1);
cscal_(&j, &signbc, &a_ref(1, j), &c__1);
/* ----------------- Begin Timing Code --------------------- */
opst += (real) ((j - 1) * 12);
/* ------------------ End Timing Code ---------------------- */
} else {
i__2 = a_subscr(j, j);
i__3 = a_subscr(j, j);
q__1.r = a[i__3].r * signbc.r - a[i__3].i * signbc.i, q__1.i =
a[i__3].r * signbc.i + a[i__3].i * signbc.r;
a[i__2].r = q__1.r, a[i__2].i = q__1.i;
}
if (ilz) {
cscal_(n, &signbc, &z___ref(1, j), &c__1);
}
/* ------------------- Begin Timing Code ---------------------- */
opst += (real) (nz * 6 + 13);
/* -------------------- End Timing Code ----------------------- */
} else {
i__2 = b_subscr(j, j);
b[i__2].r = 0.f, b[i__2].i = 0.f;
}
i__2 = j;
i__3 = a_subscr(j, j);
alpha[i__2].r = a[i__3].r, alpha[i__2].i = a[i__3].i;
i__2 = j;
i__3 = b_subscr(j, j);
beta[i__2].r = b[i__3].r, beta[i__2].i = b[i__3].i;
/* L200: */
}
/* Normal Termination */
*info = 0;
/* Exit (other than argument error) -- return optimal workspace size */
L210:
/* ---------------------- Begin Timing Code ------------------------- */
latime_1.ops += opst;
opst = 0.f;
latime_1.itcnt = (real) jiter;
/* ----------------------- End Timing Code -------------------------- */
q__1.r = (real) (*n), q__1.i = 0.f;
work[1].r = q__1.r, work[1].i = q__1.i;
return 0;
/* End of CHGEQZ */
} /* chgeqz_ */
#undef z___ref
#undef z___subscr
#undef q_ref
#undef q_subscr
#undef b_ref
#undef b_subscr
#undef a_ref
#undef a_subscr
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -