dhgeqz.c
来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 1,728 行 · 第 1/5 页
C
1,728 行
/*< A11 = A( ILAST-1, ILAST-1 ) >*/
a11 = a[ilast - 1 + (ilast - 1) * a_dim1];
/*< A21 = A( ILAST, ILAST-1 ) >*/
a21 = a[ilast + (ilast - 1) * a_dim1];
/*< A12 = A( ILAST-1, ILAST ) >*/
a12 = a[ilast - 1 + ilast * a_dim1];
/*< A22 = A( ILAST, ILAST ) >*/
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 >*/
c11r = s1 * a11 - wr * b11;
/*< C11I = -WI*B11 >*/
c11i = -wi * b11;
/*< C12 = S1*A12 >*/
c12 = s1 * a12;
/*< C21 = S1*A21 >*/
c21 = s1 * a21;
/*< C22R = S1*A22 - WR*B22 >*/
c22r = s1 * a22 - wr * b22;
/*< C22I = -WI*B22 >*/
c22i = -wi * b22;
/*< >*/
if (abs(c11r) + abs(c11i) + abs(c12) > abs(c21) + abs(c22r) + abs(
c22i)) {
/*< T = DLAPY3( C12, C11R, C11I ) >*/
t = dlapy3_(&c12, &c11r, &c11i);
/*< CZ = C12 / T >*/
cz = c12 / t;
/*< SZR = -C11R / T >*/
szr = -c11r / t;
/*< SZI = -C11I / T >*/
szi = -c11i / t;
/*< ELSE >*/
} else {
/*< CZ = DLAPY2( C22R, C22I ) >*/
cz = dlapy2_(&c22r, &c22i);
/*< IF( CZ.LE.SAFMIN ) THEN >*/
if (cz <= safmin) {
/*< CZ = ZERO >*/
cz = 0.;
/*< SZR = ONE >*/
szr = 1.;
/*< SZI = ZERO >*/
szi = 0.;
/*< ELSE >*/
} else {
/*< TEMPR = C22R / CZ >*/
tempr = c22r / cz;
/*< TEMPI = C22I / CZ >*/
tempi = c22i / cz;
/*< T = DLAPY2( CZ, C21 ) >*/
t = dlapy2_(&cz, &c21);
/*< CZ = CZ / T >*/
cz /= t;
/*< SZR = -C21*TEMPR / T >*/
szr = -c21 * tempr / t;
/*< SZI = C21*TEMPI / T >*/
szi = c21 * tempi / t;
/*< END IF >*/
}
/*< END IF >*/
}
/* Compute Givens rotation on left */
/* ( CQ SQ ) */
/* ( __ ) A or B */
/* ( -SQ CQ ) */
/*< AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 ) >*/
an = abs(a11) + abs(a12) + abs(a21) + abs(a22);
/*< BN = ABS( B11 ) + ABS( B22 ) >*/
bn = abs(b11) + abs(b22);
/*< WABS = ABS( WR ) + ABS( WI ) >*/
wabs = abs(wr) + abs(wi);
/*< IF( S1*AN.GT.WABS*BN ) THEN >*/
if (s1 * an > wabs * bn) {
/*< CQ = CZ*B11 >*/
cq = cz * b11;
/*< SQR = SZR*B22 >*/
sqr = szr * b22;
/*< SQI = -SZI*B22 >*/
sqi = -szi * b22;
/*< ELSE >*/
} else {
/*< A1R = CZ*A11 + SZR*A12 >*/
a1r = cz * a11 + szr * a12;
/*< A1I = SZI*A12 >*/
a1i = szi * a12;
/*< A2R = CZ*A21 + SZR*A22 >*/
a2r = cz * a21 + szr * a22;
/*< A2I = SZI*A22 >*/
a2i = szi * a22;
/*< CQ = DLAPY2( A1R, A1I ) >*/
cq = dlapy2_(&a1r, &a1i);
/*< IF( CQ.LE.SAFMIN ) THEN >*/
if (cq <= safmin) {
/*< CQ = ZERO >*/
cq = 0.;
/*< SQR = ONE >*/
sqr = 1.;
/*< SQI = ZERO >*/
sqi = 0.;
/*< ELSE >*/
} else {
/*< TEMPR = A1R / CQ >*/
tempr = a1r / cq;
/*< TEMPI = A1I / CQ >*/
tempi = a1i / cq;
/*< SQR = TEMPR*A2R + TEMPI*A2I >*/
sqr = tempr * a2r + tempi * a2i;
/*< SQI = TEMPI*A2R - TEMPR*A2I >*/
sqi = tempi * a2r - tempr * a2i;
/*< END IF >*/
}
/*< END IF >*/
}
/*< T = DLAPY3( CQ, SQR, SQI ) >*/
t = dlapy3_(&cq, &sqr, &sqi);
/*< CQ = CQ / T >*/
cq /= t;
/*< SQR = SQR / T >*/
sqr /= t;
/*< SQI = SQI / T >*/
sqi /= t;
/* Compute diagonal elements of QBZ */
/*< TEMPR = SQR*SZR - SQI*SZI >*/
tempr = sqr * szr - sqi * szi;
/*< TEMPI = SQR*SZI + SQI*SZR >*/
tempi = sqr * szi + sqi * szr;
/*< B1R = CQ*CZ*B11 + TEMPR*B22 >*/
b1r = cq * cz * b11 + tempr * b22;
/*< B1I = TEMPI*B22 >*/
b1i = tempi * b22;
/*< B1A = DLAPY2( B1R, B1I ) >*/
b1a = dlapy2_(&b1r, &b1i);
/*< B2R = CQ*CZ*B22 + TEMPR*B11 >*/
b2r = cq * cz * b22 + tempr * b11;
/*< B2I = -TEMPI*B11 >*/
b2i = -tempi * b11;
/*< B2A = DLAPY2( B2R, B2I ) >*/
b2a = dlapy2_(&b2r, &b2i);
/* Normalize so beta > 0, and Im( alpha1 ) > 0 */
/*< BETA( ILAST-1 ) = B1A >*/
beta[ilast - 1] = b1a;
/*< BETA( ILAST ) = B2A >*/
beta[ilast] = b2a;
/*< ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV >*/
alphar[ilast - 1] = wr * b1a * s1inv;
/*< ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV >*/
alphai[ilast - 1] = wi * b1a * s1inv;
/*< ALPHAR( ILAST ) = ( WR*B2A )*S1INV >*/
alphar[ilast] = wr * b2a * s1inv;
/*< ALPHAI( ILAST ) = -( WI*B2A )*S1INV >*/
alphai[ilast] = -(wi * b2a) * s1inv;
/* Step 3: Go to next block -- exit if finished. */
/*< ILAST = IFIRST - 1 >*/
ilast = ifirst - 1;
/*< >*/
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;
/*< ELSE >*/
} else {
/* Usual case: 3x3 or larger block, using Francis implicit */
/* double-shift */
/* 2 */
/* Eigenvalue equation is w - c w + d = 0, */
/* -1 2 -1 */
/* so compute 1st column of (A B ) - c A B + d */
/* using the formula in QZIT (from EISPACK) */
/* We assume that the block is at least 3x3 */
/*< >*/
ad11 = ascale * a[ilast - 1 + (ilast - 1) * a_dim1] / (bscale * b[
ilast - 1 + (ilast - 1) * b_dim1]);
/*< >*/
ad21 = ascale * a[ilast + (ilast - 1) * a_dim1] / (bscale * b[
ilast - 1 + (ilast - 1) * b_dim1]);
/*< >*/
ad12 = ascale * a[ilast - 1 + ilast * a_dim1] / (bscale * b[ilast
+ ilast * b_dim1]);
/*< >*/
ad22 = ascale * a[ilast + ilast * a_dim1] / (bscale * b[ilast +
ilast * b_dim1]);
/*< U12 = B( ILAST-1, ILAST ) / B( ILAST, ILAST ) >*/
u12 = b[ilast - 1 + ilast * b_dim1] / b[ilast + ilast * b_dim1];
/*< >*/
ad11l = ascale * a[ifirst + ifirst * a_dim1] / (bscale * b[ifirst
+ ifirst * b_dim1]);
/*< >*/
ad21l = ascale * a[ifirst + 1 + ifirst * a_dim1] / (bscale * b[
ifirst + ifirst * b_dim1]);
/*< >*/
ad12l = ascale * a[ifirst + (ifirst + 1) * a_dim1] / (bscale * b[
ifirst + 1 + (ifirst + 1) * b_dim1]);
/*< >*/
ad22l = ascale * a[ifirst + 1 + (ifirst + 1) * a_dim1] / (bscale *
b[ifirst + 1 + (ifirst + 1) * b_dim1]);
/*< >*/
ad32l = ascale * a[ifirst + 2 + (ifirst + 1) * a_dim1] / (bscale *
b[ifirst + 1 + (ifirst + 1) * b_dim1]);
/*< U12L = B( IFIRST, IFIRST+1 ) / B( IFIRST+1, IFIRST+1 ) >*/
u12l = b[ifirst + (ifirst + 1) * b_dim1] / b[ifirst + 1 + (ifirst
+ 1) * b_dim1];
/*< >*/
v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12
* ad11l + (ad12l - ad11l * u12l) * ad21l;
/*< >*/
v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 -
ad11l) + ad21 * u12) * ad21l;
/*< V( 3 ) = AD32L*AD21L >*/
v[2] = ad32l * ad21l;
/*< ISTART = IFIRST >*/
istart = ifirst;
/*< CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU ) >*/
dlarfg_(&c__3, v, &v[1], &c__1, &tau);
/*< V( 1 ) = ONE >*/
v[0] = 1.;
/* Sweep */
/*< DO 290 J = ISTART, ILAST - 2 >*/
i__2 = ilast - 2;
for (j = istart; j <= i__2; ++j) {
/* All but last elements: use 3x3 Householder transforms. */
/* Zero (j-1)st column of A */
/*< IF( J.GT.ISTART ) THEN >*/
if (j > istart) {
/*< V( 1 ) = A( J, J-1 ) >*/
v[0] = a[j + (j - 1) * a_dim1];
/*< V( 2 ) = A( J+1, J-1 ) >*/
v[1] = a[j + 1 + (j - 1) * a_dim1];
/*< V( 3 ) = A( J+2, J-1 ) >*/
v[2] = a[j + 2 + (j - 1) * a_dim1];
/*< CALL DLARFG( 3, A( J, J-1 ), V( 2 ), 1, TAU ) >*/
dlarfg_(&c__3, &a[j + (j - 1) * a_dim1], &v[1], &c__1, &
tau);
/*< V( 1 ) = ONE >*/
v[0] = 1.;
/*< A( J+1, J-1 ) = ZERO >*/
a[j + 1 + (j - 1) * a_dim1] = 0.;
/*< A( J+2, J-1 ) = ZERO >*/
a[j + 2 + (j - 1) * a_dim1] = 0.;
/*< END IF >*/
}
/*< DO 230 JC = J, ILASTM >*/
i__3 = ilastm;
for (jc = j; jc <= i__3; ++jc) {
/*< >*/
temp = tau * (a[j + jc * a_dim1] + v[1] * a[j + 1 + jc *
a_dim1] + v[2] * a[j + 2 + jc * a_dim1]);
/*< A( J, JC ) = A( J, JC ) - TEMP >*/
a[j + jc * a_dim1] -= temp;
/*< A( J+1, JC ) = A( J+1, JC ) - TEMP*V( 2 ) >*/
a[j + 1 + jc * a_dim1] -= temp * v[1];
/*< A( J+2, JC ) = A( J+2, JC ) - TEMP*V( 3 ) >*/
a[j + 2 + jc * a_dim1] -= temp * v[2];
/*< >*/
temp2 = tau * (b[j + jc * b_dim1] + v[1] * b[j + 1 + jc *
b_dim1] + v[2] * b[j + 2 + jc * b_dim1]);
/*< B( J, JC ) = B( J, JC ) - TEMP2 >*/
b[j + jc * b_dim1] -= temp2;
/*< B( J+1, JC ) = B( J+1, JC ) - TEMP2*V( 2 ) >*/
b[j + 1 + jc * b_dim1] -= temp2 * v[1];
/*< B( J+2, JC ) = B( J+2, JC ) - TEMP2*V( 3 ) >*/
b[j + 2 + jc * b_dim1] -= temp2 * v[2];
/*< 230 CONTINUE >*/
/* L230: */
}
/*< IF( ILQ ) THEN >*/
if (ilq) {
/*< DO 240 JR = 1, N >*/
i__3 = *n;
for (jr = 1; jr <= i__3; ++jr) {
/*< >*/
temp = tau * (q[jr + j * q_dim1] + v[1] * q[jr + (j +
1) * q_dim1] + v[2] * q[jr + (j + 2) * q_dim1]
);
/*< Q( JR, J ) = Q( JR, J ) - TEMP >*/
q[jr + j * q_dim1] -= temp;
/*< Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 ) >*/
q[jr + (j + 1) * q_dim1] -= temp * v[1];
/*< Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 ) >*/
q[jr + (j + 2) * q_dim1] -= temp * v[2];
/*< 240 CONTINUE >*/
/* L240: */
}
/*< END IF >*/
}
/* Zero j-th column of B (see DLAGBC for details) */
/* Swap rows to pivot */
/*< ILPIVT = .FALSE. >*/
ilpivt = FALSE_;
/*< TEMP = MAX( ABS( B( J+1, J+1 ) ), ABS( B( J+1, J+2 ) ) ) >*/
/* Computing MAX */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?