📄 dtgex2.c
字号:
/*< >*/
if (linfo != 0) {
goto L70;
}
/*< >*/
dorm2r_("L", "T", &m, &m, &m, tcpy, &c__4, taul, scpy, &c__4, &work[1]
, info, (ftnlen)1, (ftnlen)1);
/*< >*/
dorm2r_("R", "N", &m, &m, &m, tcpy, &c__4, taul, licop, &c__4, &work[
1], info, (ftnlen)1, (ftnlen)1);
/*< >*/
if (linfo != 0) {
goto L70;
}
/* Compute F-norm(S21) in BQRA21. (T21 is 0.) */
/*< DSCALE = ZERO >*/
dscale = 0.;
/*< DSUM = ONE >*/
dsum = 1.;
/*< DO 40 I = 1, N2 >*/
i__1 = *n2;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< CALL DLASSQ( N1, SCPY( N2+1, I ), 1, DSCALE, DSUM ) >*/
dlassq_(n1, &scpy[*n2 + 1 + (i__ << 2) - 5], &c__1, &dscale, &
dsum);
/*< 40 CONTINUE >*/
/* L40: */
}
/*< BQRA21 = DSCALE*SQRT( DSUM ) >*/
bqra21 = dscale * sqrt(dsum);
/* Decide which method to use. */
/* Weak stability test: */
/* F-norm(S21) <= O(EPS * F-norm((S, T))) */
/*< IF( BQRA21.LE.BRQA21 .AND. BQRA21.LE.THRESH ) THEN >*/
if (bqra21 <= brqa21 && bqra21 <= thresh) {
/*< CALL DLACPY( 'F', M, M, SCPY, LDST, S, LDST ) >*/
dlacpy_("F", &m, &m, scpy, &c__4, s, &c__4, (ftnlen)1);
/*< CALL DLACPY( 'F', M, M, TCPY, LDST, T, LDST ) >*/
dlacpy_("F", &m, &m, tcpy, &c__4, t, &c__4, (ftnlen)1);
/*< CALL DLACPY( 'F', M, M, IRCOP, LDST, IR, LDST ) >*/
dlacpy_("F", &m, &m, ircop, &c__4, ir, &c__4, (ftnlen)1);
/*< CALL DLACPY( 'F', M, M, LICOP, LDST, LI, LDST ) >*/
dlacpy_("F", &m, &m, licop, &c__4, li, &c__4, (ftnlen)1);
/*< ELSE IF( BRQA21.GE.THRESH ) THEN >*/
} else if (brqa21 >= thresh) {
/*< GO TO 70 >*/
goto L70;
/*< END IF >*/
}
/* Set lower triangle of B-part to zero */
/*< DO 50 I = 2, M >*/
i__1 = m;
for (i__ = 2; i__ <= i__1; ++i__) {
/*< CALL DCOPY( M-I+1, ZERO, 0, T( I, I-1 ), 1 ) >*/
i__2 = m - i__ + 1;
dcopy_(&i__2, &c_b3, &c__0, &t[i__ + ((i__ - 1) << 2) - 5], &c__1);
/*< 50 CONTINUE >*/
/* L50: */
}
/*< IF( WANDS ) THEN >*/
if (TRUE_) {
/* Strong stability test: */
/* F-norm((A-QL*S*QR', B-QL*T*QR')) <= O(EPS*F-norm((A,B))) */
/*< >*/
dlacpy_("Full", &m, &m, &a[*j1 + *j1 * a_dim1], lda, &work[m * m
+ 1], &m, (ftnlen)4);
/*< >*/
dgemm_("N", "N", &m, &m, &m, &c_b38, li, &c__4, s, &c__4, &c_b3, &
work[1], &m, (ftnlen)1, (ftnlen)1);
/*< >*/
dgemm_("N", "N", &m, &m, &m, &c_b44, &work[1], &m, ir, &c__4, &
c_b38, &work[m * m + 1], &m, (ftnlen)1, (ftnlen)1);
/*< DSCALE = ZERO >*/
dscale = 0.;
/*< DSUM = ONE >*/
dsum = 1.;
/*< CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) >*/
i__1 = m * m;
dlassq_(&i__1, &work[m * m + 1], &c__1, &dscale, &dsum);
/*< >*/
dlacpy_("Full", &m, &m, &b[*j1 + *j1 * b_dim1], ldb, &work[m * m
+ 1], &m, (ftnlen)4);
/*< >*/
dgemm_("N", "N", &m, &m, &m, &c_b38, li, &c__4, t, &c__4, &c_b3, &
work[1], &m, (ftnlen)1, (ftnlen)1);
/*< >*/
dgemm_("N", "N", &m, &m, &m, &c_b44, &work[1], &m, ir, &c__4, &
c_b38, &work[m * m + 1], &m, (ftnlen)1, (ftnlen)1);
/*< CALL DLASSQ( M*M, WORK( M*M+1 ), 1, DSCALE, DSUM ) >*/
i__1 = m * m;
dlassq_(&i__1, &work[m * m + 1], &c__1, &dscale, &dsum);
/*< SS = DSCALE*SQRT( DSUM ) >*/
ss = dscale * sqrt(dsum);
/*< DTRONG = ( SS.LE.THRESH ) >*/
dtrong = ss <= thresh;
/*< >*/
if (! dtrong) {
goto L70;
}
/*< END IF >*/
}
/* If the swap is accepted ("weakly" and "strongly"), apply the */
/* transformations and set N1-by-N2 (2,1)-block to zero. */
/*< DO 60 I = 1, N2 >*/
i__1 = *n2;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< CALL DCOPY( N1, ZERO, 0, S( N2+1, I ), 1 ) >*/
dcopy_(n1, &c_b3, &c__0, &s[*n2 + 1 + (i__ << 2) - 5], &c__1);
/*< 60 CONTINUE >*/
/* L60: */
}
/* copy back M-by-M diagonal block starting at index J1 of (A, B) */
/*< CALL DLACPY( 'F', M, M, S, LDST, A( J1, J1 ), LDA ) >*/
dlacpy_("F", &m, &m, s, &c__4, &a[*j1 + *j1 * a_dim1], lda, (ftnlen)1)
;
/*< CALL DLACPY( 'F', M, M, T, LDST, B( J1, J1 ), LDB ) >*/
dlacpy_("F", &m, &m, t, &c__4, &b[*j1 + *j1 * b_dim1], ldb, (ftnlen)1)
;
/*< CALL DCOPY( LDST*LDST, ZERO, 0, T, 1 ) >*/
dcopy_(&c__16, &c_b3, &c__0, t, &c__1);
/* Standardize existing 2-by-2 blocks. */
/*< CALL DCOPY( M*M, ZERO, 0, WORK, 1 ) >*/
i__1 = m * m;
dcopy_(&i__1, &c_b3, &c__0, &work[1], &c__1);
/*< WORK( 1 ) = ONE >*/
work[1] = 1.;
/*< T( 1, 1 ) = ONE >*/
t[0] = 1.;
/*< IDUM = LWORK - M*M - 2 >*/
idum = *lwork - m * m - 2;
/*< IF( N2.GT.1 ) THEN >*/
if (*n2 > 1) {
/*< >*/
dlagv2_(&a[*j1 + *j1 * a_dim1], lda, &b[*j1 + *j1 * b_dim1], ldb,
ar, ai, be, &work[1], &work[2], t, &t[1]);
/*< WORK( M+1 ) = -WORK( 2 ) >*/
work[m + 1] = -work[2];
/*< WORK( M+2 ) = WORK( 1 ) >*/
work[m + 2] = work[1];
/*< T( N2, N2 ) = T( 1, 1 ) >*/
t[*n2 + (*n2 << 2) - 5] = t[0];
/*< T( 1, 2 ) = -T( 2, 1 ) >*/
t[4] = -t[1];
/*< END IF >*/
}
/*< WORK( M*M ) = ONE >*/
work[m * m] = 1.;
/*< T( M, M ) = ONE >*/
t[m + (m << 2) - 5] = 1.;
/*< IF( N1.GT.1 ) THEN >*/
if (*n1 > 1) {
/*< >*/
dlagv2_(&a[*j1 + *n2 + (*j1 + *n2) * a_dim1], lda, &b[*j1 + *n2 +
(*j1 + *n2) * b_dim1], ldb, taur, taul, &work[m * m + 1],
&work[*n2 * m + *n2 + 1], &work[*n2 * m + *n2 + 2], &t[*
n2 + 1 + ((*n2 + 1) << 2) - 5], &t[m + ((m - 1) << 2) - 5]);
/*< WORK( M*M ) = WORK( N2*M+N2+1 ) >*/
work[m * m] = work[*n2 * m + *n2 + 1];
/*< WORK( M*M-1 ) = -WORK( N2*M+N2+2 ) >*/
work[m * m - 1] = -work[*n2 * m + *n2 + 2];
/*< T( M, M ) = T( N2+1, N2+1 ) >*/
t[m + (m << 2) - 5] = t[*n2 + 1 + ((*n2 + 1) << 2) - 5];
/*< T( M-1, M ) = -T( M, M-1 ) >*/
t[m - 1 + (m << 2) - 5] = -t[m + ((m - 1) << 2) - 5];
/*< END IF >*/
}
/*< >*/
dgemm_("T", "N", n2, n1, n2, &c_b38, &work[1], &m, &a[*j1 + (*j1 + *
n2) * a_dim1], lda, &c_b3, &work[m * m + 1], n2, (ftnlen)1, (
ftnlen)1);
/*< >*/
dlacpy_("Full", n2, n1, &work[m * m + 1], n2, &a[*j1 + (*j1 + *n2) *
a_dim1], lda, (ftnlen)4);
/*< >*/
dgemm_("T", "N", n2, n1, n2, &c_b38, &work[1], &m, &b[*j1 + (*j1 + *
n2) * b_dim1], ldb, &c_b3, &work[m * m + 1], n2, (ftnlen)1, (
ftnlen)1);
/*< >*/
dlacpy_("Full", n2, n1, &work[m * m + 1], n2, &b[*j1 + (*j1 + *n2) *
b_dim1], ldb, (ftnlen)4);
/*< >*/
dgemm_("N", "N", &m, &m, &m, &c_b38, li, &c__4, &work[1], &m, &c_b3, &
work[m * m + 1], &m, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', M, M, WORK( M*M+1 ), M, LI, LDST ) >*/
dlacpy_("Full", &m, &m, &work[m * m + 1], &m, li, &c__4, (ftnlen)4);
/*< >*/
dgemm_("N", "N", n2, n1, n1, &c_b38, &a[*j1 + (*j1 + *n2) * a_dim1],
lda, &t[*n2 + 1 + ((*n2 + 1) << 2) - 5], &c__4, &c_b3, &work[1],
n2, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', N2, N1, WORK, N2, A( J1, J1+N2 ), LDA ) >*/
dlacpy_("Full", n2, n1, &work[1], n2, &a[*j1 + (*j1 + *n2) * a_dim1],
lda, (ftnlen)4);
/*< >*/
dgemm_("N", "N", n2, n1, n1, &c_b38, &b[*j1 + (*j1 + *n2) * b_dim1],
lda, &t[*n2 + 1 + ((*n2 + 1) << 2) - 5], &c__4, &c_b3, &work[1],
n2, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', N2, N1, WORK, N2, B( J1, J1+N2 ), LDB ) >*/
dlacpy_("Full", n2, n1, &work[1], n2, &b[*j1 + (*j1 + *n2) * b_dim1],
ldb, (ftnlen)4);
/*< >*/
dgemm_("T", "N", &m, &m, &m, &c_b38, ir, &c__4, t, &c__4, &c_b3, &
work[1], &m, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', M, M, WORK, M, IR, LDST ) >*/
dlacpy_("Full", &m, &m, &work[1], &m, ir, &c__4, (ftnlen)4);
/* Accumulate transformations into Q and Z if requested. */
/*< IF( WANTQ ) THEN >*/
if (*wantq) {
/*< >*/
dgemm_("N", "N", n, &m, &m, &c_b38, &q[*j1 * q_dim1 + 1], ldq, li,
&c__4, &c_b3, &work[1], n, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', N, M, WORK, N, Q( 1, J1 ), LDQ ) >*/
dlacpy_("Full", n, &m, &work[1], n, &q[*j1 * q_dim1 + 1], ldq, (
ftnlen)4);
/*< END IF >*/
}
/*< IF( WANTZ ) THEN >*/
if (*wantz) {
/*< >*/
dgemm_("N", "N", n, &m, &m, &c_b38, &z__[*j1 * z_dim1 + 1], ldz,
ir, &c__4, &c_b3, &work[1], n, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', N, M, WORK, N, Z( 1, J1 ), LDZ ) >*/
dlacpy_("Full", n, &m, &work[1], n, &z__[*j1 * z_dim1 + 1], ldz, (
ftnlen)4);
/*< END IF >*/
}
/* Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and */
/* (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)). */
/*< I = J1 + M >*/
i__ = *j1 + m;
/*< IF( I.LE.N ) THEN >*/
if (i__ <= *n) {
/*< >*/
i__1 = *n - i__ + 1;
dgemm_("T", "N", &m, &i__1, &m, &c_b38, li, &c__4, &a[*j1 + i__ *
a_dim1], lda, &c_b3, &work[1], &m, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', M, N-I+1, WORK, M, A( J1, I ), LDA ) >*/
i__1 = *n - i__ + 1;
dlacpy_("Full", &m, &i__1, &work[1], &m, &a[*j1 + i__ * a_dim1],
lda, (ftnlen)4);
/*< >*/
i__1 = *n - i__ + 1;
dgemm_("T", "N", &m, &i__1, &m, &c_b38, li, &c__4, &b[*j1 + i__ *
b_dim1], lda, &c_b3, &work[1], &m, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', M, N-I+1, WORK, M, B( J1, I ), LDA ) >*/
i__1 = *n - i__ + 1;
dlacpy_("Full", &m, &i__1, &work[1], &m, &b[*j1 + i__ * b_dim1],
lda, (ftnlen)4);
/*< END IF >*/
}
/*< I = J1 - 1 >*/
i__ = *j1 - 1;
/*< IF( I.GT.0 ) THEN >*/
if (i__ > 0) {
/*< >*/
dgemm_("N", "N", &i__, &m, &m, &c_b38, &a[*j1 * a_dim1 + 1], lda,
ir, &c__4, &c_b3, &work[1], &i__, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', I, M, WORK, I, A( 1, J1 ), LDA ) >*/
dlacpy_("Full", &i__, &m, &work[1], &i__, &a[*j1 * a_dim1 + 1],
lda, (ftnlen)4);
/*< >*/
dgemm_("N", "N", &i__, &m, &m, &c_b38, &b[*j1 * b_dim1 + 1], ldb,
ir, &c__4, &c_b3, &work[1], &i__, (ftnlen)1, (ftnlen)1);
/*< CALL DLACPY( 'Full', I, M, WORK, I, B( 1, J1 ), LDB ) >*/
dlacpy_("Full", &i__, &m, &work[1], &i__, &b[*j1 * b_dim1 + 1],
ldb, (ftnlen)4);
/*< END IF >*/
}
/* Exit with INFO = 0 if swap was successfully performed. */
/*< RETURN >*/
return 0;
/*< END IF >*/
}
/* Exit with INFO = 1 if swap was rejected. */
/*< 70 CONTINUE >*/
L70:
/*< INFO = 1 >*/
*info = 1;
/*< RETURN >*/
return 0;
/* End of DTGEX2 */
/*< END >*/
} /* dtgex2_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -