📄 zhgeqz.c
字号:
ad11.r = z__1.r, ad11.i = z__1.i;
i__2 = a_subscr(ilast, ilast - 1);
z__2.r = ascale * a[i__2].r, z__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast - 1, ilast - 1);
z__3.r = bscale * b[i__3].r, z__3.i = bscale * b[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad21.r = z__1.r, ad21.i = z__1.i;
i__2 = a_subscr(ilast - 1, ilast);
z__2.r = ascale * a[i__2].r, z__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast, ilast);
z__3.r = bscale * b[i__3].r, z__3.i = bscale * b[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad12.r = z__1.r, ad12.i = z__1.i;
i__2 = a_subscr(ilast, ilast);
z__2.r = ascale * a[i__2].r, z__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast, ilast);
z__3.r = bscale * b[i__3].r, z__3.i = bscale * b[i__3].i;
z_div(&z__1, &z__2, &z__3);
ad22.r = z__1.r, ad22.i = z__1.i;
z__2.r = u12.r * ad21.r - u12.i * ad21.i, z__2.i = u12.r * ad21.i
+ u12.i * ad21.r;
z__1.r = ad22.r - z__2.r, z__1.i = ad22.i - z__2.i;
abi22.r = z__1.r, abi22.i = z__1.i;
z__2.r = ad11.r + abi22.r, z__2.i = ad11.i + abi22.i;
z__1.r = z__2.r * .5, z__1.i = z__2.i * .5;
t.r = z__1.r, t.i = z__1.i;
pow_zi(&z__4, &t, &c__2);
z__5.r = ad12.r * ad21.r - ad12.i * ad21.i, z__5.i = ad12.r *
ad21.i + ad12.i * ad21.r;
z__3.r = z__4.r + z__5.r, z__3.i = z__4.i + z__5.i;
z__6.r = ad11.r * ad22.r - ad11.i * ad22.i, z__6.i = ad11.r *
ad22.i + ad11.i * ad22.r;
z__2.r = z__3.r - z__6.r, z__2.i = z__3.i - z__6.i;
z_sqrt(&z__1, &z__2);
rtdisc.r = z__1.r, rtdisc.i = z__1.i;
z__1.r = t.r - abi22.r, z__1.i = t.i - abi22.i;
z__2.r = t.r - abi22.r, z__2.i = t.i - abi22.i;
temp = z__1.r * rtdisc.r + d_imag(&z__2) * d_imag(&rtdisc);
if (temp <= 0.) {
z__1.r = t.r + rtdisc.r, z__1.i = t.i + rtdisc.i;
shift.r = z__1.r, shift.i = z__1.i;
} else {
z__1.r = t.r - rtdisc.r, z__1.i = t.i - rtdisc.i;
shift.r = z__1.r, shift.i = z__1.i;
}
/* ------------------- Begin Timing Code ---------------------- */
opst += 116.;
/* -------------------- End Timing Code ----------------------- */
} else {
/* Exceptional shift. Chosen for no particularly good reason. */
i__2 = a_subscr(ilast - 1, ilast);
z__4.r = ascale * a[i__2].r, z__4.i = ascale * a[i__2].i;
i__3 = b_subscr(ilast - 1, ilast - 1);
z__5.r = bscale * b[i__3].r, z__5.i = bscale * b[i__3].i;
z_div(&z__3, &z__4, &z__5);
d_cnjg(&z__2, &z__3);
z__1.r = eshift.r + z__2.r, z__1.i = eshift.i + z__2.i;
eshift.r = z__1.r, eshift.i = z__1.i;
shift.r = eshift.r, shift.i = eshift.i;
/* ------------------- Begin Timing Code ---------------------- */
opst += 15.;
/* -------------------- 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);
z__2.r = ascale * a[i__3].r, z__2.i = ascale * a[i__3].i;
i__4 = b_subscr(j, j);
z__4.r = bscale * b[i__4].r, z__4.i = bscale * b[i__4].i;
z__3.r = shift.r * z__4.r - shift.i * z__4.i, z__3.i = shift.r *
z__4.i + shift.i * z__4.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
temp = (d__1 = ctemp.r, abs(d__1)) + (d__2 = d_imag(&ctemp), abs(
d__2));
i__3 = a_subscr(j + 1, j);
temp2 = ascale * ((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&
a_ref(j + 1, j)), abs(d__2)));
tempr = max(temp,temp2);
if (tempr < 1. && tempr != 0.) {
temp /= tempr;
temp2 /= tempr;
}
i__3 = a_subscr(j, j - 1);
if (((d__1 = a[i__3].r, abs(d__1)) + (d__2 = d_imag(&a_ref(j, j -
1)), abs(d__2))) * temp2 <= temp * atol) {
goto L90;
}
/* L80: */
}
istart = ifirst;
i__2 = a_subscr(ifirst, ifirst);
z__2.r = ascale * a[i__2].r, z__2.i = ascale * a[i__2].i;
i__3 = b_subscr(ifirst, ifirst);
z__4.r = bscale * b[i__3].r, z__4.i = bscale * b[i__3].i;
z__3.r = shift.r * z__4.r - shift.i * z__4.i, z__3.i = shift.r *
z__4.i + shift.i * z__4.r;
z__1.r = z__2.r - z__3.r, z__1.i = z__2.i - z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
/* --------------------- Begin Timing Code ----------------------- */
opst += -6.;
/* ---------------------- End Timing Code ------------------------ */
L90:
/* Do an implicit-shift QZ sweep.
Initial Q */
i__2 = a_subscr(istart + 1, istart);
z__1.r = ascale * a[i__2].r, z__1.i = ascale * a[i__2].i;
ctemp2.r = z__1.r, ctemp2.i = z__1.i;
/* --------------------- Begin Timing Code ----------------------- */
opst += (doublereal) ((ilast - istart) * 18 + 2);
/* ---------------------- End Timing Code ------------------------ */
zlartg_(&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;
zlartg_(&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., a[i__3].i = 0.;
}
i__3 = ilastm;
for (jc = j; jc <= i__3; ++jc) {
i__4 = a_subscr(j, jc);
z__2.r = c__ * a[i__4].r, z__2.i = c__ * a[i__4].i;
i__5 = a_subscr(j + 1, jc);
z__3.r = s.r * a[i__5].r - s.i * a[i__5].i, z__3.i = s.r * a[
i__5].i + s.i * a[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = a_subscr(j + 1, jc);
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = a_subscr(j, jc);
z__2.r = z__3.r * a[i__5].r - z__3.i * a[i__5].i, z__2.i =
z__3.r * a[i__5].i + z__3.i * a[i__5].r;
i__6 = a_subscr(j + 1, jc);
z__5.r = c__ * a[i__6].r, z__5.i = c__ * a[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
a[i__4].r = z__1.r, a[i__4].i = z__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);
z__2.r = c__ * b[i__4].r, z__2.i = c__ * b[i__4].i;
i__5 = b_subscr(j + 1, jc);
z__3.r = s.r * b[i__5].r - s.i * b[i__5].i, z__3.i = s.r * b[
i__5].i + s.i * b[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp2.r = z__1.r, ctemp2.i = z__1.i;
i__4 = b_subscr(j + 1, jc);
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = b_subscr(j, jc);
z__2.r = z__3.r * b[i__5].r - z__3.i * b[i__5].i, z__2.i =
z__3.r * b[i__5].i + z__3.i * b[i__5].r;
i__6 = b_subscr(j + 1, jc);
z__5.r = c__ * b[i__6].r, z__5.i = c__ * b[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
b[i__4].r = z__1.r, b[i__4].i = z__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);
z__2.r = c__ * q[i__4].r, z__2.i = c__ * q[i__4].i;
d_cnjg(&z__4, &s);
i__5 = q_subscr(jr, j + 1);
z__3.r = z__4.r * q[i__5].r - z__4.i * q[i__5].i, z__3.i =
z__4.r * q[i__5].i + z__4.i * q[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = q_subscr(jr, j + 1);
z__3.r = -s.r, z__3.i = -s.i;
i__5 = q_subscr(jr, j);
z__2.r = z__3.r * q[i__5].r - z__3.i * q[i__5].i, z__2.i =
z__3.r * q[i__5].i + z__3.i * q[i__5].r;
i__6 = q_subscr(jr, j + 1);
z__4.r = c__ * q[i__6].r, z__4.i = c__ * q[i__6].i;
z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
q[i__4].r = z__1.r, q[i__4].i = z__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;
zlartg_(&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., b[i__3].i = 0.;
/* 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);
z__2.r = c__ * a[i__4].r, z__2.i = c__ * a[i__4].i;
i__5 = a_subscr(jr, j);
z__3.r = s.r * a[i__5].r - s.i * a[i__5].i, z__3.i = s.r * a[
i__5].i + s.i * a[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = a_subscr(jr, j);
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = a_subscr(jr, j + 1);
z__2.r = z__3.r * a[i__5].r - z__3.i * a[i__5].i, z__2.i =
z__3.r * a[i__5].i + z__3.i * a[i__5].r;
i__6 = a_subscr(jr, j);
z__5.r = c__ * a[i__6].r, z__5.i = c__ * a[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
a[i__4].r = z__1.r, a[i__4].i = z__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);
z__2.r = c__ * b[i__4].r, z__2.i = c__ * b[i__4].i;
i__5 = b_subscr(jr, j);
z__3.r = s.r * b[i__5].r - s.i * b[i__5].i, z__3.i = s.r * b[
i__5].i + s.i * b[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = b_subscr(jr, j);
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = b_subscr(jr, j + 1);
z__2.r = z__3.r * b[i__5].r - z__3.i * b[i__5].i, z__2.i =
z__3.r * b[i__5].i + z__3.i * b[i__5].r;
i__6 = b_subscr(jr, j);
z__5.r = c__ * b[i__6].r, z__5.i = c__ * b[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
b[i__4].r = z__1.r, b[i__4].i = z__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);
z__2.r = c__ * z__[i__4].r, z__2.i = c__ * z__[i__4].i;
i__5 = z___subscr(jr, j);
z__3.r = s.r * z__[i__5].r - s.i * z__[i__5].i, z__3.i =
s.r * z__[i__5].i + s.i * z__[i__5].r;
z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__4 = z___subscr(jr, j);
d_cnjg(&z__4, &s);
z__3.r = -z__4.r, z__3.i = -z__4.i;
i__5 = z___subscr(jr, j + 1);
z__2.r = z__3.r * z__[i__5].r - z__3.i * z__[i__5].i,
z__2.i = z__3.r * z__[i__5].i + z__3.i * z__[i__5]
.r;
i__6 = z___subscr(jr, j);
z__5.r = c__ * z__[i__6].r, z__5.i = c__ * z__[i__6].i;
z__1.r = z__2.r + z__5.r, z__1.i = z__2.i + z__5.i;
z__[i__4].r = z__1.r, z__[i__4].i = z__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 += (doublereal) (ilast - istart) * (doublereal) ((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.;
/* ---------------------- End Timing Code ------------------------
L170: */
}
/* Drop-through = non-convergence */
L180:
*info = ilast;
/* ---------------------- Begin Timing Code ------------------------- */
latime_1.ops += opst;
opst = 0.;
/* ----------------------- 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 = z_abs(&b_ref(j, j));
if (absb > safmin) {
i__2 = b_subscr(j, j);
z__2.r = b[i__2].r / absb, z__2.i = b[i__2].i / absb;
d_cnjg(&z__1, &z__2);
signbc.r = z__1.r, signbc.i = z__1.i;
i__2 = b_subscr(j, j);
b[i__2].r = absb, b[i__2].i = 0.;
if (ilschr) {
i__2 = j - 1;
zscal_(&i__2, &signbc, &b_ref(1, j), &c__1);
zscal_(&j, &signbc, &a_ref(1, j), &c__1);
/* ----------------- Begin Timing Code --------------------- */
opst += (doublereal) ((j - 1) * 12);
/* ------------------ End Timing Code ---------------------- */
} else {
i__2 = a_subscr(j, j);
i__3 = a_subscr(j, j);
z__1.r = a[i__3].r * signbc.r - a[i__3].i * signbc.i, z__1.i =
a[i__3].r * signbc.i + a[i__3].i * signbc.r;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
}
if (ilz) {
zscal_(n, &signbc, &z___ref(1, j), &c__1);
}
/* ------------------- Begin Timing Code ---------------------- */
opst += (doublereal) (nz * 6 + 13);
/* -------------------- End Timing Code ----------------------- */
} else {
i__2 = b_subscr(j, j);
b[i__2].r = 0., b[i__2].i = 0.;
}
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.;
latime_1.itcnt = (doublereal) jiter;
/* ----------------------- End Timing Code -------------------------- */
z__1.r = (doublereal) (*n), z__1.i = 0.;
work[1].r = z__1.r, work[1].i = z__1.i;
return 0;
/* End of ZHGEQZ */
} /* zhgeqz_ */
#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 + -