📄 dhgeqz.c
字号:
}
/* Reset counters */
iiter = 0;
eshift = 0.;
if (! ilschr) {
ilastm = ilast;
if (ifrstm > ilast) {
ifrstm = *ilo;
}
}
goto L350;
/* QZ step */
/* This iteration only involves rows/columns IFIRST:ILAST. We */
/* assume IFIRST < ILAST, and that the diagonal of B is non-zero. */
L110:
++iiter;
if (! ilschr) {
ifrstm = ifirst;
}
/* Compute single shifts. */
/* At this point, IFIRST < ILAST, and the diagonal elements of */
/* B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in */
/* magnitude) */
if (iiter / 10 * 10 == iiter) {
/* Exceptional shift. Chosen for no particularly good reason. */
/* (Single shift only.) */
if ((doublereal) maxit * safmin * abs(a[ilast - 1 + ilast * a_dim1])
< abs(b[ilast - 1 + (ilast - 1) * b_dim1])) {
eshift += a[ilast - 1 + ilast * a_dim1] / b[ilast - 1 + (ilast - 1) * b_dim1];
} else {
eshift += 1. / (safmin * (doublereal) maxit);
}
s1 = 1.;
wr = eshift;
} else {
/* Shifts based on the generalized eigenvalues of the */
/* bottom-right 2x2 block of A and B. The first eigenvalue */
/* returned by DLAG2 is the Wilkinson shift (AEP p.512), */
d__1 = safmin * 100.;
dlag2_(&a[ilast - 1 + (ilast - 1) * a_dim1], lda, &b[ilast - 1 + (ilast - 1) * b_dim1],
ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);
temp = max(s1, safmin * max(max(1.,abs(wr)),abs(wi)));
if (wi != 0.) {
goto L200;
}
}
/* Fiddle with shift to avoid overflow */
temp = min(ascale,1.) * (safmax * .5);
if (s1 > temp) {
scale = temp / s1;
} else {
scale = 1.;
}
temp = min(bscale,1.) * (safmax * .5);
if (abs(wr) > temp) {
scale = min(scale, temp/abs(wr));
}
s1 *= scale;
wr *= scale;
/* Now check for two consecutive small subdiagonals. */
for (j = ilast - 1; j > ifirst; --j) {
istart = j;
temp = abs(s1 * a[j + (j - 1) * a_dim1]);
temp2 = abs(s1 * a[j + j * a_dim1] - wr * b[j + j * b_dim1]);
tempr = max(temp,temp2);
if (tempr < 1. && tempr != 0.) {
temp /= tempr;
temp2 /= tempr;
}
if (abs(ascale * a[j + 1 + j * a_dim1] * temp) <= ascale * atol * temp2) {
goto L130;
}
}
istart = ifirst;
L130:
/* Do an implicit single-shift QZ sweep. */
/* Initial Q */
temp = s1 * a[istart + istart * a_dim1] - wr * b[istart + istart * b_dim1];
temp2 = s1 * a[istart + 1 + istart * a_dim1];
dlartg_(&temp, &temp2, &c, &s, &tempr);
/* Sweep */
for (j = istart; j < ilast; ++j) {
if (j > istart) {
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;
}
}
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 <= min(j+2, 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 <= j; ++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;
}
}
}
goto L350;
/* Use Francis double-shift */
/* Note: the Francis double-shift should work with real shifts, */
/* but only if the block is at least 3x3. */
/* This code may break if this point is reached with */
/* a 2x2 block with real eigenvalues. */
L200:
if (ifirst + 1 == ilast) {
/* Special case -- 2x2 block with complex eigenvectors */
/* Step 1: Standardize, that is, rotate so that */
/* ( B11 0 ) */
/* B = ( ) with B11 non-negative. */
/* ( 0 B22 ) */
dlasv2_(&b[ilast - 1 + (ilast - 1) * b_dim1], &b[ilast - 1 + ilast * b_dim1],
&b[ilast + ilast * b_dim1], &b22, &b11, &sr, &cr, &sl, &cl);
if (b11 < 0.) {
cr = -cr;
sr = -sr;
b11 = -b11;
b22 = -b22;
}
i__1 = ilastm + 1 - ifirst;
drot_(&i__1, &a[ilast - 1 + (ilast - 1) * a_dim1], lda, &a[ilast + (ilast - 1) * a_dim1], lda, &cl, &sl);
i__1 = ilast + 1 - ifrstm;
drot_(&i__1, &a[ifrstm + (ilast - 1) * a_dim1], &c__1, &a[ifrstm + ilast * a_dim1], &c__1, &cr, &sr);
if (ilast < ilastm) {
i__1 = ilastm - ilast;
drot_(&i__1, &b[ilast - 1 + (ilast + 1) * b_dim1], ldb, &b[ilast + (ilast + 1) * b_dim1], lda, &cl, &sl);
}
if (ifrstm < ilast - 1) {
i__1 = ifirst - ifrstm;
drot_(&i__1, &b[ifrstm + (ilast - 1) * b_dim1], &c__1, &b[ifrstm + ilast * b_dim1], &c__1, &cr, &sr);
}
if (ilq) {
drot_(n, &q[(ilast - 1) * q_dim1 + 1], &c__1, &q[ilast * q_dim1 + 1], &c__1, &cl, &sl);
}
if (ilz) {
drot_(n, &z[(ilast - 1) * z_dim1 + 1], &c__1, &z[ilast * z_dim1 + 1], &c__1, &cr, &sr);
}
b[ilast - 1 + (ilast - 1) * b_dim1] = b11;
b[ilast - 1 + ilast * b_dim1] = 0.;
b[ilast + (ilast - 1) * b_dim1] = 0.;
b[ilast + ilast * b_dim1] = b22;
/* If B22 is negative, negate column ILAST */
if (b22 < 0.) {
for (j = ifrstm; j <= ilast; ++j) {
a[j + ilast * a_dim1] = -a[j + ilast * a_dim1];
b[j + ilast * b_dim1] = -b[j + ilast * b_dim1];
}
if (ilz) {
for (j = 1; j <= *n; ++j) {
z[j + ilast * z_dim1] = -z[j + ilast * z_dim1];
}
}
}
/* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) */
/* Recompute shift */
d__1 = safmin * 100.;
dlag2_(&a[ilast - 1 + (ilast - 1) * a_dim1], lda, &b[ilast - 1 + (ilast - 1) * b_dim1],
ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);
/* If standardization has perturbed the shift onto real line, */
/* do another (real single-shift) QR step. */
if (wi == 0.) {
goto L350;
}
s1inv = 1. / s1;
/* Do EISPACK (QZVAL) computation of alpha and beta */
a11 = a[ilast - 1 + (ilast - 1) * a_dim1];
a21 = a[ilast + (ilast - 1) * a_dim1];
a12 = a[ilast - 1 + ilast * a_dim1];
a22 = a[ilast + ilast * a_dim1];
/* Compute complex Givens rotation on right */
/* (Assume some element of C = (sA - wB) > unfl ) */
/* __ */
/* (sA - wB) ( CZ -SZ ) */
/* ( SZ CZ ) */
c11r = s1 * a11 - wr * b11;
c11i = -wi * b11;
c12 = s1 * a12;
c21 = s1 * a21;
c22r = s1 * a22 - wr * b22;
c22i = -wi * b22;
if (abs(c11r) + abs(c11i) + abs(c12) > abs(c21) + abs(c22r) + abs(c22i)) {
t = dlapy3_(&c12, &c11r, &c11i);
cz = c12 / t;
szr = -c11r / t;
szi = -c11i / t;
} else {
cz = dlapy2_(&c22r, &c22i);
if (cz <= safmin) {
cz = 0.;
szr = 1.;
szi = 0.;
} else {
tempr = c22r / cz;
tempi = c22i / cz;
t = dlapy2_(&cz, &c21);
cz /= t;
szr = -c21 * tempr / t;
szi = c21 * tempi / t;
}
}
/* Compute Givens rotation on left */
/* ( CQ SQ ) */
/* ( __ ) A or B */
/* ( -SQ CQ ) */
an = abs(a11) + abs(a12) + abs(a21) + abs(a22);
bn = abs(b11) + abs(b22);
wabs = abs(wr) + abs(wi);
if (s1 * an > wabs * bn) {
cq = cz * b11;
sqr = szr * b22;
sqi = -szi * b22;
} else {
a1r = cz * a11 + szr * a12;
a1i = szi * a12;
a2r = cz * a21 + szr * a22;
a2i = szi * a22;
cq = dlapy2_(&a1r, &a1i);
if (cq <= safmin) {
cq = 0.;
sqr = 1.;
sqi = 0.;
} else {
tempr = a1r / cq;
tempi = a1i / cq;
sqr = tempr * a2r + tempi * a2i;
sqi = tempi * a2r - tempr * a2i;
}
}
t = dlapy3_(&cq, &sqr, &sqi);
cq /= t;
sqr /= t;
sqi /= t;
/* Compute diagonal elements of QBZ */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -