dhgeqz.c
来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 1,728 行 · 第 1/5 页
C
1,728 行
/*< >*/
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( ONE, ABS( WR ), ABS( WI ) ) ) >*/
/* Computing MAX */
/* Computing MAX */
d__3 = 1., d__4 = abs(wr), d__3 = max(d__3,d__4), d__4 = abs(wi);
d__1 = s1, d__2 = safmin * max(d__3,d__4);
temp = max(d__1,d__2);
/*< >*/
if (wi != 0.) {
goto L200;
}
/*< END IF >*/
}
/* Fiddle with shift to avoid overflow */
/*< TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX ) >*/
temp = min(ascale,1.) * (safmax * .5);
/*< IF( S1.GT.TEMP ) THEN >*/
if (s1 > temp) {
/*< SCALE = TEMP / S1 >*/
scale = temp / s1;
/*< ELSE >*/
} else {
/*< SCALE = ONE >*/
scale = 1.;
/*< END IF >*/
}
/*< TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX ) >*/
temp = min(bscale,1.) * (safmax * .5);
/*< >*/
if (abs(wr) > temp) {
/* Computing MIN */
d__1 = scale, d__2 = temp / abs(wr);
scale = min(d__1,d__2);
}
/*< S1 = SCALE*S1 >*/
s1 = scale * s1;
/*< WR = SCALE*WR >*/
wr = scale * wr;
/* Now check for two consecutive small subdiagonals. */
/*< DO 120 J = ILAST - 1, IFIRST + 1, -1 >*/
i__2 = ifirst + 1;
for (j = ilast - 1; j >= i__2; --j) {
/*< ISTART = J >*/
istart = j;
/*< TEMP = ABS( S1*A( J, J-1 ) ) >*/
temp = (d__1 = s1 * a[j + (j - 1) * a_dim1], abs(d__1));
/*< TEMP2 = ABS( S1*A( J, J )-WR*B( J, J ) ) >*/
temp2 = (d__1 = s1 * a[j + j * a_dim1] - wr * b[j + j * b_dim1],
abs(d__1));
/*< TEMPR = MAX( TEMP, TEMP2 ) >*/
tempr = max(temp,temp2);
/*< IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN >*/
if (tempr < 1. && tempr != 0.) {
/*< TEMP = TEMP / TEMPR >*/
temp /= tempr;
/*< TEMP2 = TEMP2 / TEMPR >*/
temp2 /= tempr;
/*< END IF >*/
}
/*< >*/
if ((d__1 = ascale * a[j + 1 + j * a_dim1] * temp, abs(d__1)) <=
ascale * atol * temp2) {
goto L130;
}
/*< 120 CONTINUE >*/
/* L120: */
}
/*< ISTART = IFIRST >*/
istart = ifirst;
/*< 130 CONTINUE >*/
L130:
/* Do an implicit single-shift QZ sweep. */
/* Initial Q */
/*< TEMP = S1*A( ISTART, ISTART ) - WR*B( ISTART, ISTART ) >*/
temp = s1 * a[istart + istart * a_dim1] - wr * b[istart + istart *
b_dim1];
/*< TEMP2 = S1*A( ISTART+1, ISTART ) >*/
temp2 = s1 * a[istart + 1 + istart * a_dim1];
/*< CALL DLARTG( TEMP, TEMP2, C, S, TEMPR ) >*/
dlartg_(&temp, &temp2, &c__, &s, &tempr);
/* Sweep */
/*< DO 190 J = ISTART, ILAST - 1 >*/
i__2 = ilast - 1;
for (j = istart; j <= i__2; ++j) {
/*< IF( J.GT.ISTART ) THEN >*/
if (j > istart) {
/*< TEMP = A( J, J-1 ) >*/
temp = a[j + (j - 1) * a_dim1];
/*< CALL DLARTG( TEMP, A( J+1, J-1 ), C, S, A( J, J-1 ) ) >*/
dlartg_(&temp, &a[j + 1 + (j - 1) * a_dim1], &c__, &s, &a[j +
(j - 1) * a_dim1]);
/*< A( J+1, J-1 ) = ZERO >*/
a[j + 1 + (j - 1) * a_dim1] = 0.;
/*< END IF >*/
}
/*< DO 140 JC = J, ILASTM >*/
i__3 = ilastm;
for (jc = j; jc <= i__3; ++jc) {
/*< TEMP = C*A( J, JC ) + S*A( J+1, JC ) >*/
temp = c__ * a[j + jc * a_dim1] + s * a[j + 1 + jc * a_dim1];
/*< A( J+1, JC ) = -S*A( J, JC ) + C*A( J+1, JC ) >*/
a[j + 1 + jc * a_dim1] = -s * a[j + jc * a_dim1] + c__ * a[j
+ 1 + jc * a_dim1];
/*< A( J, JC ) = TEMP >*/
a[j + jc * a_dim1] = temp;
/*< TEMP2 = C*B( J, JC ) + S*B( J+1, JC ) >*/
temp2 = c__ * b[j + jc * b_dim1] + s * b[j + 1 + jc * b_dim1];
/*< B( J+1, JC ) = -S*B( J, JC ) + C*B( J+1, JC ) >*/
b[j + 1 + jc * b_dim1] = -s * b[j + jc * b_dim1] + c__ * b[j
+ 1 + jc * b_dim1];
/*< B( J, JC ) = TEMP2 >*/
b[j + jc * b_dim1] = temp2;
/*< 140 CONTINUE >*/
/* L140: */
}
/*< IF( ILQ ) THEN >*/
if (ilq) {
/*< DO 150 JR = 1, N >*/
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
/*< TEMP = C*Q( JR, J ) + S*Q( JR, J+1 ) >*/
temp = c__ * q[jr + j * q_dim1] + s * q[jr + (j + 1) *
q_dim1];
/*< Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 ) >*/
q[jr + (j + 1) * q_dim1] = -s * q[jr + j * q_dim1] + c__ *
q[jr + (j + 1) * q_dim1];
/*< Q( JR, J ) = TEMP >*/
q[jr + j * q_dim1] = temp;
/*< 150 CONTINUE >*/
/* L150: */
}
/*< END IF >*/
}
/*< TEMP = B( J+1, J+1 ) >*/
temp = b[j + 1 + (j + 1) * b_dim1];
/*< CALL DLARTG( TEMP, B( J+1, J ), C, S, B( J+1, J+1 ) ) >*/
dlartg_(&temp, &b[j + 1 + j * b_dim1], &c__, &s, &b[j + 1 + (j +
1) * b_dim1]);
/*< B( J+1, J ) = ZERO >*/
b[j + 1 + j * b_dim1] = 0.;
/*< DO 160 JR = IFRSTM, MIN( J+2, ILAST ) >*/
/* Computing MIN */
i__4 = j + 2;
i__3 = min(i__4,ilast);
for (jr = ifrstm; jr <= i__3; ++jr) {
/*< TEMP = C*A( JR, J+1 ) + S*A( JR, J ) >*/
temp = c__ * a[jr + (j + 1) * a_dim1] + s * a[jr + j * a_dim1]
;
/*< A( JR, J ) = -S*A( JR, J+1 ) + C*A( JR, J ) >*/
a[jr + j * a_dim1] = -s * a[jr + (j + 1) * a_dim1] + c__ * a[
jr + j * a_dim1];
/*< A( JR, J+1 ) = TEMP >*/
a[jr + (j + 1) * a_dim1] = temp;
/*< 160 CONTINUE >*/
/* L160: */
}
/*< DO 170 JR = IFRSTM, J >*/
i__3 = j;
for (jr = ifrstm; jr <= i__3; ++jr) {
/*< TEMP = C*B( JR, J+1 ) + S*B( JR, J ) >*/
temp = c__ * b[jr + (j + 1) * b_dim1] + s * b[jr + j * b_dim1]
;
/*< B( JR, J ) = -S*B( JR, J+1 ) + C*B( JR, J ) >*/
b[jr + j * b_dim1] = -s * b[jr + (j + 1) * b_dim1] + c__ * b[
jr + j * b_dim1];
/*< B( JR, J+1 ) = TEMP >*/
b[jr + (j + 1) * b_dim1] = temp;
/*< 170 CONTINUE >*/
/* L170: */
}
/*< IF( ILZ ) THEN >*/
if (ilz) {
/*< DO 180 JR = 1, N >*/
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
/*< TEMP = C*Z( JR, J+1 ) + S*Z( JR, J ) >*/
temp = c__ * z__[jr + (j + 1) * z_dim1] + s * z__[jr + j *
z_dim1];
/*< Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J ) >*/
z__[jr + j * z_dim1] = -s * z__[jr + (j + 1) * z_dim1] +
c__ * z__[jr + j * z_dim1];
/*< Z( JR, J+1 ) = TEMP >*/
z__[jr + (j + 1) * z_dim1] = temp;
/*< 180 CONTINUE >*/
/* L180: */
}
/*< END IF >*/
}
/*< 190 CONTINUE >*/
/* L190: */
}
/*< GO TO 350 >*/
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. */
/*< 200 CONTINUE >*/
L200:
/*< IF( IFIRST+1.EQ.ILAST ) THEN >*/
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.LT.ZERO ) THEN >*/
if (b11 < 0.) {
/*< CR = -CR >*/
cr = -cr;
/*< SR = -SR >*/
sr = -sr;
/*< B11 = -B11 >*/
b11 = -b11;
/*< B22 = -B22 >*/
b22 = -b22;
/*< END IF >*/
}
/*< >*/
i__2 = ilastm + 1 - ifirst;
drot_(&i__2, &a[ilast - 1 + (ilast - 1) * a_dim1], lda, &a[ilast
+ (ilast - 1) * a_dim1], lda, &cl, &sl);
/*< >*/
i__2 = ilast + 1 - ifrstm;
drot_(&i__2, &a[ifrstm + (ilast - 1) * a_dim1], &c__1, &a[ifrstm
+ ilast * a_dim1], &c__1, &cr, &sr);
/*< >*/
if (ilast < ilastm) {
i__2 = ilastm - ilast;
drot_(&i__2, &b[ilast - 1 + (ilast + 1) * b_dim1], ldb, &b[
ilast + (ilast + 1) * b_dim1], lda, &cl, &sl);
}
/*< >*/
if (ifrstm < ilast - 1) {
i__2 = ifirst - ifrstm;
drot_(&i__2, &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 ) = B11 >*/
b[ilast - 1 + (ilast - 1) * b_dim1] = b11;
/*< B( ILAST-1, ILAST ) = ZERO >*/
b[ilast - 1 + ilast * b_dim1] = 0.;
/*< B( ILAST, ILAST-1 ) = ZERO >*/
b[ilast + (ilast - 1) * b_dim1] = 0.;
/*< B( ILAST, ILAST ) = B22 >*/
b[ilast + ilast * b_dim1] = b22;
/* If B22 is negative, negate column ILAST */
/*< IF( B22.LT.ZERO ) THEN >*/
if (b22 < 0.) {
/*< DO 210 J = IFRSTM, ILAST >*/
i__2 = ilast;
for (j = ifrstm; j <= i__2; ++j) {
/*< A( J, ILAST ) = -A( J, ILAST ) >*/
a[j + ilast * a_dim1] = -a[j + ilast * a_dim1];
/*< B( J, ILAST ) = -B( J, ILAST ) >*/
b[j + ilast * b_dim1] = -b[j + ilast * b_dim1];
/*< 210 CONTINUE >*/
/* L210: */
}
/*< IF( ILZ ) THEN >*/
if (ilz) {
/*< DO 220 J = 1, N >*/
i__2 = *n;
for (j = 1; j <= i__2; ++j) {
/*< Z( J, ILAST ) = -Z( J, ILAST ) >*/
z__[j + ilast * z_dim1] = -z__[j + ilast * z_dim1];
/*< 220 CONTINUE >*/
/* L220: */
}
/*< END IF >*/
}
/*< END IF >*/
}
/* 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 = ONE / S1 >*/
s1inv = 1. / s1;
/* Do EISPACK (QZVAL) computation of alpha and beta */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?