📄 zlahqr.c
字号:
z_sqrt(&z__1, &z__2);
y.r = z__1.r, y.i = z__1.i;
/*< >*/
if (x.r * y.r + d_imag(&x) * d_imag(&y) < 0.) {
z__1.r = -y.r, z__1.i = -y.i;
y.r = z__1.r, y.i = z__1.i;
}
/*< T = T - ZLADIV( U, ( X+Y ) ) >*/
z__3.r = x.r + y.r, z__3.i = x.i + y.i;
zladiv_(&z__2, &u, &z__3);
z__1.r = t.r - z__2.r, z__1.i = t.i - z__2.i;
t.r = z__1.r, t.i = z__1.i;
/*< END IF >*/
}
/*< END IF >*/
}
/* Look for two consecutive small subdiagonal elements. */
/*< DO 40 M = I - 1, L + 1, -1 >*/
i__2 = l + 1;
for (m = i__ - 1; m >= i__2; --m) {
/* Determine the effect of starting the single-shift QR */
/* iteration at row M, and see if this would make H(M,M-1) */
/* negligible. */
/*< H11 = H( M, M ) >*/
i__3 = m + m * h_dim1;
h11.r = h__[i__3].r, h11.i = h__[i__3].i;
/*< H22 = H( M+1, M+1 ) >*/
i__3 = m + 1 + (m + 1) * h_dim1;
h22.r = h__[i__3].r, h22.i = h__[i__3].i;
/*< H11S = H11 - T >*/
z__1.r = h11.r - t.r, z__1.i = h11.i - t.i;
h11s.r = z__1.r, h11s.i = z__1.i;
/*< H21 = H( M+1, M ) >*/
i__3 = m + 1 + m * h_dim1;
h21 = h__[i__3].r;
/*< S = CABS1( H11S ) + ABS( H21 ) >*/
s = (d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(d__2))
+ abs(h21);
/*< H11S = H11S / S >*/
z__1.r = h11s.r / s, z__1.i = h11s.i / s;
h11s.r = z__1.r, h11s.i = z__1.i;
/*< H21 = H21 / S >*/
h21 /= s;
/*< V( 1 ) = H11S >*/
v[0].r = h11s.r, v[0].i = h11s.i;
/*< V( 2 ) = H21 >*/
v[1].r = h21, v[1].i = 0.;
/*< H10 = H( M, M-1 ) >*/
i__3 = m + (m - 1) * h_dim1;
h10 = h__[i__3].r;
/*< TST1 = CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) ) >*/
tst1 = ((d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(
d__2))) * ((d__3 = h11.r, abs(d__3)) + (d__4 = d_imag(&
h11), abs(d__4)) + ((d__5 = h22.r, abs(d__5)) + (d__6 =
d_imag(&h22), abs(d__6))));
/*< >*/
if ((d__1 = h10 * h21, abs(d__1)) <= ulp * tst1) {
goto L50;
}
/*< 40 CONTINUE >*/
/* L40: */
}
/*< H11 = H( L, L ) >*/
i__2 = l + l * h_dim1;
h11.r = h__[i__2].r, h11.i = h__[i__2].i;
/*< H22 = H( L+1, L+1 ) >*/
i__2 = l + 1 + (l + 1) * h_dim1;
h22.r = h__[i__2].r, h22.i = h__[i__2].i;
/*< H11S = H11 - T >*/
z__1.r = h11.r - t.r, z__1.i = h11.i - t.i;
h11s.r = z__1.r, h11s.i = z__1.i;
/*< H21 = H( L+1, L ) >*/
i__2 = l + 1 + l * h_dim1;
h21 = h__[i__2].r;
/*< S = CABS1( H11S ) + ABS( H21 ) >*/
s = (d__1 = h11s.r, abs(d__1)) + (d__2 = d_imag(&h11s), abs(d__2)) +
abs(h21);
/*< H11S = H11S / S >*/
z__1.r = h11s.r / s, z__1.i = h11s.i / s;
h11s.r = z__1.r, h11s.i = z__1.i;
/*< H21 = H21 / S >*/
h21 /= s;
/*< V( 1 ) = H11S >*/
v[0].r = h11s.r, v[0].i = h11s.i;
/*< V( 2 ) = H21 >*/
v[1].r = h21, v[1].i = 0.;
/*< 50 CONTINUE >*/
L50:
/* Single-shift QR step */
/*< DO 100 K = M, I - 1 >*/
i__2 = i__ - 1;
for (k = m; k <= i__2; ++k) {
/* The first iteration of this loop determines a reflection G */
/* from the vector V and applies it from left and right to H, */
/* thus creating a nonzero bulge below the subdiagonal. */
/* Each subsequent iteration determines a reflection G to */
/* restore the Hessenberg form in the (K-1)th column, and thus */
/* chases the bulge one step toward the bottom of the active */
/* submatrix. */
/* V(2) is always real before the call to ZLARFG, and hence */
/* after the call T2 ( = T1*V(2) ) is also real. */
/*< >*/
if (k > m) {
zcopy_(&c__2, &h__[k + (k - 1) * h_dim1], &c__1, v, &c__1);
}
/*< CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 ) >*/
zlarfg_(&c__2, v, &v[1], &c__1, &t1);
/*< IF( K.GT.M ) THEN >*/
if (k > m) {
/*< H( K, K-1 ) = V( 1 ) >*/
i__3 = k + (k - 1) * h_dim1;
h__[i__3].r = v[0].r, h__[i__3].i = v[0].i;
/*< H( K+1, K-1 ) = ZERO >*/
i__3 = k + 1 + (k - 1) * h_dim1;
h__[i__3].r = 0., h__[i__3].i = 0.;
/*< END IF >*/
}
/*< V2 = V( 2 ) >*/
v2.r = v[1].r, v2.i = v[1].i;
/*< T2 = DBLE( T1*V2 ) >*/
z__1.r = t1.r * v2.r - t1.i * v2.i, z__1.i = t1.r * v2.i + t1.i *
v2.r;
t2 = z__1.r;
/* Apply G from the left to transform the rows of the matrix */
/* in columns K to I2. */
/*< DO 60 J = K, I2 >*/
i__3 = i2;
for (j = k; j <= i__3; ++j) {
/*< SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J ) >*/
d_cnjg(&z__3, &t1);
i__4 = k + j * h_dim1;
z__2.r = z__3.r * h__[i__4].r - z__3.i * h__[i__4].i, z__2.i =
z__3.r * h__[i__4].i + z__3.i * h__[i__4].r;
i__5 = k + 1 + j * h_dim1;
z__4.r = t2 * h__[i__5].r, z__4.i = t2 * h__[i__5].i;
z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
sum.r = z__1.r, sum.i = z__1.i;
/*< H( K, J ) = H( K, J ) - SUM >*/
i__4 = k + j * h_dim1;
i__5 = k + j * h_dim1;
z__1.r = h__[i__5].r - sum.r, z__1.i = h__[i__5].i - sum.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/*< H( K+1, J ) = H( K+1, J ) - SUM*V2 >*/
i__4 = k + 1 + j * h_dim1;
i__5 = k + 1 + j * h_dim1;
z__2.r = sum.r * v2.r - sum.i * v2.i, z__2.i = sum.r * v2.i +
sum.i * v2.r;
z__1.r = h__[i__5].r - z__2.r, z__1.i = h__[i__5].i - z__2.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/*< 60 CONTINUE >*/
/* L60: */
}
/* Apply G from the right to transform the columns of the */
/* matrix in rows I1 to min(K+2,I). */
/*< DO 70 J = I1, MIN( K+2, I ) >*/
/* Computing MIN */
i__4 = k + 2;
i__3 = min(i__4,i__);
for (j = i1; j <= i__3; ++j) {
/*< SUM = T1*H( J, K ) + T2*H( J, K+1 ) >*/
i__4 = j + k * h_dim1;
z__2.r = t1.r * h__[i__4].r - t1.i * h__[i__4].i, z__2.i =
t1.r * h__[i__4].i + t1.i * h__[i__4].r;
i__5 = j + (k + 1) * h_dim1;
z__3.r = t2 * h__[i__5].r, z__3.i = t2 * h__[i__5].i;
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;
/*< H( J, K ) = H( J, K ) - SUM >*/
i__4 = j + k * h_dim1;
i__5 = j + k * h_dim1;
z__1.r = h__[i__5].r - sum.r, z__1.i = h__[i__5].i - sum.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/*< H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 ) >*/
i__4 = j + (k + 1) * h_dim1;
i__5 = j + (k + 1) * h_dim1;
d_cnjg(&z__3, &v2);
z__2.r = sum.r * z__3.r - sum.i * z__3.i, z__2.i = sum.r *
z__3.i + sum.i * z__3.r;
z__1.r = h__[i__5].r - z__2.r, z__1.i = h__[i__5].i - z__2.i;
h__[i__4].r = z__1.r, h__[i__4].i = z__1.i;
/*< 70 CONTINUE >*/
/* L70: */
}
/*< IF( WANTZ ) THEN >*/
if (*wantz) {
/* Accumulate transformations in the matrix Z */
/*< DO 80 J = ILOZ, IHIZ >*/
i__3 = *ihiz;
for (j = *iloz; j <= i__3; ++j) {
/*< SUM = T1*Z( J, K ) + T2*Z( J, K+1 ) >*/
i__4 = j + k * z_dim1;
z__2.r = t1.r * z__[i__4].r - t1.i * z__[i__4].i, z__2.i =
t1.r * z__[i__4].i + t1.i * z__[i__4].r;
i__5 = j + (k + 1) * z_dim1;
z__3.r = t2 * z__[i__5].r, z__3.i = t2 * z__[i__5].i;
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;
/*< Z( J, K ) = Z( J, K ) - SUM >*/
i__4 = j + k * z_dim1;
i__5 = j + k * z_dim1;
z__1.r = z__[i__5].r - sum.r, z__1.i = z__[i__5].i -
sum.i;
z__[i__4].r = z__1.r, z__[i__4].i = z__1.i;
/*< Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 ) >*/
i__4 = j + (k + 1) * z_dim1;
i__5 = j + (k + 1) * z_dim1;
d_cnjg(&z__3, &v2);
z__2.r = sum.r * z__3.r - sum.i * z__3.i, z__2.i = sum.r *
z__3.i + sum.i * z__3.r;
z__1.r = z__[i__5].r - z__2.r, z__1.i = z__[i__5].i -
z__2.i;
z__[i__4].r = z__1.r, z__[i__4].i = z__1.i;
/*< 80 CONTINUE >*/
/* L80: */
}
/*< END IF >*/
}
/*< IF( K.EQ.M .AND. M.GT.L ) THEN >*/
if (k == m && m > l) {
/* If the QR step was started at row M > L because two */
/* consecutive small subdiagonals were found, then extra */
/* scaling must be performed to ensure that H(M,M-1) remains */
/* real. */
/*< TEMP = ONE - T1 >*/
z__1.r = 1. - t1.r, z__1.i = 0. - t1.i;
temp.r = z__1.r, temp.i = z__1.i;
/*< TEMP = TEMP / ABS( TEMP ) >*/
d__1 = z_abs(&temp);
z__1.r = temp.r / d__1, z__1.i = temp.i / d__1;
temp.r = z__1.r, temp.i = z__1.i;
/*< H( M+1, M ) = H( M+1, M )*DCONJG( TEMP ) >*/
i__3 = m + 1 + m * h_dim1;
i__4 = m + 1 + m * h_dim1;
d_cnjg(&z__2, &temp);
z__1.r = h__[i__4].r * z__2.r - h__[i__4].i * z__2.i, z__1.i =
h__[i__4].r * z__2.i + h__[i__4].i * z__2.r;
h__[i__3].r = z__1.r, h__[i__3].i = z__1.i;
/*< >*/
if (m + 2 <= i__) {
i__3 = m + 2 + (m + 1) * h_dim1;
i__4 = m + 2 + (m + 1) * h_dim1;
z__1.r = h__[i__4].r * temp.r - h__[i__4].i * temp.i,
z__1.i = h__[i__4].r * temp.i + h__[i__4].i *
temp.r;
h__[i__3].r = z__1.r, h__[i__3].i = z__1.i;
}
/*< DO 90 J = M, I >*/
i__3 = i__;
for (j = m; j <= i__3; ++j) {
/*< IF( J.NE.M+1 ) THEN >*/
if (j != m + 1) {
/*< >*/
if (i2 > j) {
i__4 = i2 - j;
zscal_(&i__4, &temp, &h__[j + (j + 1) * h_dim1],
ldh);
}
/*< CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 ) >*/
i__4 = j - i1;
d_cnjg(&z__1, &temp);
zscal_(&i__4, &z__1, &h__[i1 + j * h_dim1], &c__1);
/*< IF( WANTZ ) THEN >*/
if (*wantz) {
/*< >*/
d_cnjg(&z__1, &temp);
zscal_(&nz, &z__1, &z__[*iloz + j * z_dim1], &
c__1);
/*< END IF >*/
}
/*< END IF >*/
}
/*< 90 CONTINUE >*/
/* L90: */
}
/*< END IF >*/
}
/*< 100 CONTINUE >*/
/* L100: */
}
/* Ensure that H(I,I-1) is real. */
/*< TEMP = H( I, I-1 ) >*/
i__2 = i__ + (i__ - 1) * h_dim1;
temp.r = h__[i__2].r, temp.i = h__[i__2].i;
/*< IF( DIMAG( TEMP ).NE.RZERO ) THEN >*/
if (d_imag(&temp) != 0.) {
/*< RTEMP = ABS( TEMP ) >*/
rtemp = z_abs(&temp);
/*< H( I, I-1 ) = RTEMP >*/
i__2 = i__ + (i__ - 1) * h_dim1;
h__[i__2].r = rtemp, h__[i__2].i = 0.;
/*< TEMP = TEMP / RTEMP >*/
z__1.r = temp.r / rtemp, z__1.i = temp.i / rtemp;
temp.r = z__1.r, temp.i = z__1.i;
/*< >*/
if (i2 > i__) {
i__2 = i2 - i__;
d_cnjg(&z__1, &temp);
zscal_(&i__2, &z__1, &h__[i__ + (i__ + 1) * h_dim1], ldh);
}
/*< CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 ) >*/
i__2 = i__ - i1;
zscal_(&i__2, &temp, &h__[i1 + i__ * h_dim1], &c__1);
/*< IF( WANTZ ) THEN >*/
if (*wantz) {
/*< CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 ) >*/
zscal_(&nz, &temp, &z__[*iloz + i__ * z_dim1], &c__1);
/*< END IF >*/
}
/*< END IF >*/
}
/*< 110 CONTINUE >*/
/* L110: */
}
/* Failure to converge in remaining number of iterations */
/*< INFO = I >*/
*info = i__;
/*< RETURN >*/
return 0;
/*< 120 CONTINUE >*/
L120:
/* H(I,I-1) is negligible: one eigenvalue has converged. */
/*< W( I ) = H( I, I ) >*/
i__1 = i__;
i__2 = i__ + i__ * h_dim1;
w[i__1].r = h__[i__2].r, w[i__1].i = h__[i__2].i;
/* Decrement number of remaining iterations, and return to start of */
/* the main loop with new value of I. */
/*< ITN = ITN - ITS >*/
itn -= its;
/*< I = L - 1 >*/
i__ = l - 1;
/*< GO TO 10 >*/
goto L10;
/*< 130 CONTINUE >*/
L130:
/*< RETURN >*/
return 0;
/* End of ZLAHQR */
/*< END >*/
} /* zlahqr_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -