dhgeqz.c
来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 1,728 行 · 第 1/5 页
C
1,728 行
/*< ELSE >*/
} else {
/*< ILAZRO = .FALSE. >*/
ilazro = FALSE_;
/*< END IF >*/
}
/*< END IF >*/
}
/* Test 2: for B(j,j)=0 */
/*< IF( ABS( B( J, J ) ).LT.BTOL ) THEN >*/
if ((d__1 = b[j + j * b_dim1], abs(d__1)) < btol) {
/*< B( J, J ) = ZERO >*/
b[j + j * b_dim1] = 0.;
/* Test 1a: Check for 2 consecutive small subdiagonals in A */
/*< ILAZR2 = .FALSE. >*/
ilazr2 = FALSE_;
/*< IF( .NOT.ILAZRO ) THEN >*/
if (! ilazro) {
/*< TEMP = ABS( A( J, J-1 ) ) >*/
temp = (d__1 = a[j + (j - 1) * a_dim1], abs(d__1));
/*< TEMP2 = ABS( A( J, J ) ) >*/
temp2 = (d__1 = a[j + j * a_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 (temp * (ascale * (d__1 = a[j + 1 + j * a_dim1], abs(
d__1))) <= temp2 * (ascale * atol)) {
ilazr2 = TRUE_;
}
/*< END IF >*/
}
/* If both tests pass (1 & 2), i.e., the leading diagonal */
/* element of B in the block is zero, split a 1x1 block off */
/* at the top. (I.e., at the J-th row/column) The leading */
/* diagonal element of the remainder can also be zero, so */
/* this may have to be done repeatedly. */
/*< IF( ILAZRO .OR. ILAZR2 ) THEN >*/
if (ilazro || ilazr2) {
/*< DO 40 JCH = J, ILAST - 1 >*/
i__3 = ilast - 1;
for (jch = j; jch <= i__3; ++jch) {
/*< TEMP = A( JCH, JCH ) >*/
temp = a[jch + jch * a_dim1];
/*< >*/
dlartg_(&temp, &a[jch + 1 + jch * a_dim1], &c__, &s, &
a[jch + jch * a_dim1]);
/*< A( JCH+1, JCH ) = ZERO >*/
a[jch + 1 + jch * a_dim1] = 0.;
/*< >*/
i__4 = ilastm - jch;
drot_(&i__4, &a[jch + (jch + 1) * a_dim1], lda, &a[
jch + 1 + (jch + 1) * a_dim1], lda, &c__, &s);
/*< >*/
i__4 = ilastm - jch;
drot_(&i__4, &b[jch + (jch + 1) * b_dim1], ldb, &b[
jch + 1 + (jch + 1) * b_dim1], ldb, &c__, &s);
/*< >*/
if (ilq) {
drot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1)
* q_dim1 + 1], &c__1, &c__, &s);
}
/*< >*/
if (ilazr2) {
a[jch + (jch - 1) * a_dim1] *= c__;
}
/*< ILAZR2 = .FALSE. >*/
ilazr2 = FALSE_;
/*< IF( ABS( B( JCH+1, JCH+1 ) ).GE.BTOL ) THEN >*/
if ((d__1 = b[jch + 1 + (jch + 1) * b_dim1], abs(d__1)
) >= btol) {
/*< IF( JCH+1.GE.ILAST ) THEN >*/
if (jch + 1 >= ilast) {
/*< GO TO 80 >*/
goto L80;
/*< ELSE >*/
} else {
/*< IFIRST = JCH + 1 >*/
ifirst = jch + 1;
/*< GO TO 110 >*/
goto L110;
/*< END IF >*/
}
/*< END IF >*/
}
/*< B( JCH+1, JCH+1 ) = ZERO >*/
b[jch + 1 + (jch + 1) * b_dim1] = 0.;
/*< 40 CONTINUE >*/
/* L40: */
}
/*< GO TO 70 >*/
goto L70;
/*< ELSE >*/
} else {
/* Only test 2 passed -- chase the zero to B(ILAST,ILAST) */
/* Then process as in the case B(ILAST,ILAST)=0 */
/*< DO 50 JCH = J, ILAST - 1 >*/
i__3 = ilast - 1;
for (jch = j; jch <= i__3; ++jch) {
/*< TEMP = B( JCH, JCH+1 ) >*/
temp = b[jch + (jch + 1) * b_dim1];
/*< >*/
dlartg_(&temp, &b[jch + 1 + (jch + 1) * b_dim1], &c__,
&s, &b[jch + (jch + 1) * b_dim1]);
/*< B( JCH+1, JCH+1 ) = ZERO >*/
b[jch + 1 + (jch + 1) * b_dim1] = 0.;
/*< >*/
if (jch < ilastm - 1) {
i__4 = ilastm - jch - 1;
drot_(&i__4, &b[jch + (jch + 2) * b_dim1], ldb, &
b[jch + 1 + (jch + 2) * b_dim1], ldb, &
c__, &s);
}
/*< >*/
i__4 = ilastm - jch + 2;
drot_(&i__4, &a[jch + (jch - 1) * a_dim1], lda, &a[
jch + 1 + (jch - 1) * a_dim1], lda, &c__, &s);
/*< >*/
if (ilq) {
drot_(n, &q[jch * q_dim1 + 1], &c__1, &q[(jch + 1)
* q_dim1 + 1], &c__1, &c__, &s);
}
/*< TEMP = A( JCH+1, JCH ) >*/
temp = a[jch + 1 + jch * a_dim1];
/*< >*/
dlartg_(&temp, &a[jch + 1 + (jch - 1) * a_dim1], &c__,
&s, &a[jch + 1 + jch * a_dim1]);
/*< A( JCH+1, JCH-1 ) = ZERO >*/
a[jch + 1 + (jch - 1) * a_dim1] = 0.;
/*< >*/
i__4 = jch + 1 - ifrstm;
drot_(&i__4, &a[ifrstm + jch * a_dim1], &c__1, &a[
ifrstm + (jch - 1) * a_dim1], &c__1, &c__, &s)
;
/*< >*/
i__4 = jch - ifrstm;
drot_(&i__4, &b[ifrstm + jch * b_dim1], &c__1, &b[
ifrstm + (jch - 1) * b_dim1], &c__1, &c__, &s)
;
/*< >*/
if (ilz) {
drot_(n, &z__[jch * z_dim1 + 1], &c__1, &z__[(jch
- 1) * z_dim1 + 1], &c__1, &c__, &s);
}
/*< 50 CONTINUE >*/
/* L50: */
}
/*< GO TO 70 >*/
goto L70;
/*< END IF >*/
}
/*< ELSE IF( ILAZRO ) THEN >*/
} else if (ilazro) {
/* Only test 1 passed -- work on J:ILAST */
/*< IFIRST = J >*/
ifirst = j;
/*< GO TO 110 >*/
goto L110;
/*< END IF >*/
}
/* Neither test passed -- try next J */
/*< 60 CONTINUE >*/
/* L60: */
}
/* (Drop-through is "impossible") */
/*< INFO = N + 1 >*/
*info = *n + 1;
/*< GO TO 420 >*/
goto L420;
/* B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a */
/* 1x1 block. */
/*< 70 CONTINUE >*/
L70:
/*< TEMP = A( ILAST, ILAST ) >*/
temp = a[ilast + ilast * a_dim1];
/*< >*/
dlartg_(&temp, &a[ilast + (ilast - 1) * a_dim1], &c__, &s, &a[ilast +
ilast * a_dim1]);
/*< A( ILAST, ILAST-1 ) = ZERO >*/
a[ilast + (ilast - 1) * a_dim1] = 0.;
/*< >*/
i__2 = ilast - ifrstm;
drot_(&i__2, &a[ifrstm + ilast * a_dim1], &c__1, &a[ifrstm + (ilast -
1) * a_dim1], &c__1, &c__, &s);
/*< >*/
i__2 = ilast - ifrstm;
drot_(&i__2, &b[ifrstm + ilast * b_dim1], &c__1, &b[ifrstm + (ilast -
1) * b_dim1], &c__1, &c__, &s);
/*< >*/
if (ilz) {
drot_(n, &z__[ilast * z_dim1 + 1], &c__1, &z__[(ilast - 1) *
z_dim1 + 1], &c__1, &c__, &s);
}
/* A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, */
/* and BETA */
/*< 80 CONTINUE >*/
L80:
/*< IF( B( ILAST, ILAST ).LT.ZERO ) THEN >*/
if (b[ilast + ilast * b_dim1] < 0.) {
/*< IF( ILSCHR ) THEN >*/
if (ilschr) {
/*< DO 90 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];
/*< 90 CONTINUE >*/
/* L90: */
}
/*< ELSE >*/
} else {
/*< A( ILAST, ILAST ) = -A( ILAST, ILAST ) >*/
a[ilast + ilast * a_dim1] = -a[ilast + ilast * a_dim1];
/*< B( ILAST, ILAST ) = -B( ILAST, ILAST ) >*/
b[ilast + ilast * b_dim1] = -b[ilast + ilast * b_dim1];
/*< END IF >*/
}
/*< IF( ILZ ) THEN >*/
if (ilz) {
/*< DO 100 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];
/*< 100 CONTINUE >*/
/* L100: */
}
/*< END IF >*/
}
/*< END IF >*/
}
/*< ALPHAR( ILAST ) = A( ILAST, ILAST ) >*/
alphar[ilast] = a[ilast + ilast * a_dim1];
/*< ALPHAI( ILAST ) = ZERO >*/
alphai[ilast] = 0.;
/*< BETA( ILAST ) = B( ILAST, ILAST ) >*/
beta[ilast] = b[ilast + ilast * b_dim1];
/* Go to next block -- exit if finished. */
/*< ILAST = ILAST - 1 >*/
--ilast;
/*< >*/
if (ilast < *ilo) {
goto L380;
}
/* Reset counters */
/*< IITER = 0 >*/
iiter = 0;
/*< ESHIFT = ZERO >*/
eshift = 0.;
/*< IF( .NOT.ILSCHR ) THEN >*/
if (! ilschr) {
/*< ILASTM = ILAST >*/
ilastm = ilast;
/*< >*/
if (ifrstm > ilast) {
ifrstm = *ilo;
}
/*< END IF >*/
}
/*< GO TO 350 >*/
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. */
/*< 110 CONTINUE >*/
L110:
/*< IITER = IITER + 1 >*/
++iiter;
/*< IF( .NOT.ILSCHR ) THEN >*/
if (! ilschr) {
/*< IFRSTM = IFIRST >*/
ifrstm = ifirst;
/*< END IF >*/
}
/* 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.EQ.IITER ) THEN >*/
if (iiter / 10 * 10 == iiter) {
/* Exceptional shift. Chosen for no particularly good reason. */
/* (Single shift only.) */
/*< >*/
if ((doublereal) maxit * safmin * (d__1 = a[ilast - 1 + ilast *
a_dim1], abs(d__1)) < (d__2 = b[ilast - 1 + (ilast - 1) *
b_dim1], abs(d__2))) {
/*< >*/
eshift += a[ilast - 1 + ilast * a_dim1] / b[ilast - 1 + (
ilast - 1) * b_dim1];
/*< ELSE >*/
} else {
/*< ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) ) >*/
eshift += 1. / (safmin * (doublereal) maxit);
/*< END IF >*/
}
/*< S1 = ONE >*/
s1 = 1.;
/*< WR = ESHIFT >*/
wr = eshift;
/*< ELSE >*/
} 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), */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?