📄 dtgsen.c
字号:
/*< PL = ZERO >*/
*pl = 0.;
/*< PR = ZERO >*/
*pr = 0.;
/*< END IF >*/
}
/*< IF( WANTD ) THEN >*/
if (wantd) {
/*< DIF( 1 ) = ZERO >*/
dif[1] = 0.;
/*< DIF( 2 ) = ZERO >*/
dif[2] = 0.;
/*< END IF >*/
}
/*< GO TO 60 >*/
goto L60;
/*< END IF >*/
}
/*< >*/
if (pair) {
++ks;
}
/*< END IF >*/
}
/*< END IF >*/
}
/*< 30 CONTINUE >*/
/* L30: */
}
/*< IF( WANTP ) THEN >*/
if (wantp) {
/* Solve generalized Sylvester equation for R and L */
/* and compute PL and PR. */
/*< N1 = M >*/
n1 = *m;
/*< N2 = N - M >*/
n2 = *n - *m;
/*< I = N1 + 1 >*/
i__ = n1 + 1;
/*< IJB = 0 >*/
ijb = 0;
/*< CALL DLACPY( 'Full', N1, N2, A( 1, I ), LDA, WORK, N1 ) >*/
dlacpy_("Full", &n1, &n2, &a[i__ * a_dim1 + 1], lda, &work[1], &n1, (
ftnlen)4);
/*< >*/
dlacpy_("Full", &n1, &n2, &b[i__ * b_dim1 + 1], ldb, &work[n1 * n2 +
1], &n1, (ftnlen)4);
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i__ + i__ * a_dim1]
, lda, &work[1], &n1, &b[b_offset], ldb, &b[i__ + i__ *
b_dim1], ldb, &work[n1 * n2 + 1], &n1, &dscale, &dif[1], &
work[(n1 * n2 << 1) + 1], &i__1, &iwork[1], &ierr, (ftnlen)1);
/* Estimate the reciprocal of norms of "projections" onto left */
/* and right eigenspaces. */
/*< RDSCAL = ZERO >*/
rdscal = 0.;
/*< DSUM = ONE >*/
dsum = 1.;
/*< CALL DLASSQ( N1*N2, WORK, 1, RDSCAL, DSUM ) >*/
i__1 = n1 * n2;
dlassq_(&i__1, &work[1], &c__1, &rdscal, &dsum);
/*< PL = RDSCAL*SQRT( DSUM ) >*/
*pl = rdscal * sqrt(dsum);
/*< IF( PL.EQ.ZERO ) THEN >*/
if (*pl == 0.) {
/*< PL = ONE >*/
*pl = 1.;
/*< ELSE >*/
} else {
/*< PL = DSCALE / ( SQRT( DSCALE*DSCALE / PL+PL )*SQRT( PL ) ) >*/
*pl = dscale / (sqrt(dscale * dscale / *pl + *pl) * sqrt(*pl));
/*< END IF >*/
}
/*< RDSCAL = ZERO >*/
rdscal = 0.;
/*< DSUM = ONE >*/
dsum = 1.;
/*< CALL DLASSQ( N1*N2, WORK( N1*N2+1 ), 1, RDSCAL, DSUM ) >*/
i__1 = n1 * n2;
dlassq_(&i__1, &work[n1 * n2 + 1], &c__1, &rdscal, &dsum);
/*< PR = RDSCAL*SQRT( DSUM ) >*/
*pr = rdscal * sqrt(dsum);
/*< IF( PR.EQ.ZERO ) THEN >*/
if (*pr == 0.) {
/*< PR = ONE >*/
*pr = 1.;
/*< ELSE >*/
} else {
/*< PR = DSCALE / ( SQRT( DSCALE*DSCALE / PR+PR )*SQRT( PR ) ) >*/
*pr = dscale / (sqrt(dscale * dscale / *pr + *pr) * sqrt(*pr));
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( WANTD ) THEN >*/
if (wantd) {
/* Compute estimates of Difu and Difl. */
/*< IF( WANTD1 ) THEN >*/
if (wantd1) {
/*< N1 = M >*/
n1 = *m;
/*< N2 = N - M >*/
n2 = *n - *m;
/*< I = N1 + 1 >*/
i__ = n1 + 1;
/*< IJB = IDIFJB >*/
ijb = 3;
/* Frobenius norm-based Difu-estimate. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i__ + i__ *
a_dim1], lda, &work[1], &n1, &b[b_offset], ldb, &b[i__ +
i__ * b_dim1], ldb, &work[n1 * n2 + 1], &n1, &dscale, &
dif[1], &work[(n1 << 1) * n2 + 1], &i__1, &iwork[1], &
ierr, (ftnlen)1);
/* Frobenius norm-based Difl-estimate. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("N", &ijb, &n2, &n1, &a[i__ + i__ * a_dim1], lda, &a[
a_offset], lda, &work[1], &n2, &b[i__ + i__ * b_dim1],
ldb, &b[b_offset], ldb, &work[n1 * n2 + 1], &n2, &dscale,
&dif[2], &work[(n1 << 1) * n2 + 1], &i__1, &iwork[1], &
ierr, (ftnlen)1);
/*< ELSE >*/
} else {
/* Compute 1-norm-based estimates of Difu and Difl using */
/* reversed communication with DLACON. In each step a */
/* generalized Sylvester equation or a transposed variant */
/* is solved. */
/*< KASE = 0 >*/
kase = 0;
/*< N1 = M >*/
n1 = *m;
/*< N2 = N - M >*/
n2 = *n - *m;
/*< I = N1 + 1 >*/
i__ = n1 + 1;
/*< IJB = 0 >*/
ijb = 0;
/*< MN2 = 2*N1*N2 >*/
mn2 = (n1 << 1) * n2;
/* 1-norm-based estimate of Difu. */
/*< 40 CONTINUE >*/
L40:
/*< >*/
dlacon_(&mn2, &work[mn2 + 1], &work[1], &iwork[1], &dif[1], &kase)
;
/*< IF( KASE.NE.0 ) THEN >*/
if (kase != 0) {
/*< IF( KASE.EQ.1 ) THEN >*/
if (kase == 1) {
/* Solve generalized Sylvester equation. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("N", &ijb, &n1, &n2, &a[a_offset], lda, &a[i__ +
i__ * a_dim1], lda, &work[1], &n1, &b[b_offset],
ldb, &b[i__ + i__ * b_dim1], ldb, &work[n1 * n2 +
1], &n1, &dscale, &dif[1], &work[(n1 << 1) * n2 +
1], &i__1, &iwork[1], &ierr, (ftnlen)1);
/*< ELSE >*/
} else {
/* Solve the transposed variant. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("T", &ijb, &n1, &n2, &a[a_offset], lda, &a[i__ +
i__ * a_dim1], lda, &work[1], &n1, &b[b_offset],
ldb, &b[i__ + i__ * b_dim1], ldb, &work[n1 * n2 +
1], &n1, &dscale, &dif[1], &work[(n1 << 1) * n2 +
1], &i__1, &iwork[1], &ierr, (ftnlen)1);
/*< END IF >*/
}
/*< GO TO 40 >*/
goto L40;
/*< END IF >*/
}
/*< DIF( 1 ) = DSCALE / DIF( 1 ) >*/
dif[1] = dscale / dif[1];
/* 1-norm-based estimate of Difl. */
/*< 50 CONTINUE >*/
L50:
/*< >*/
dlacon_(&mn2, &work[mn2 + 1], &work[1], &iwork[1], &dif[2], &kase)
;
/*< IF( KASE.NE.0 ) THEN >*/
if (kase != 0) {
/*< IF( KASE.EQ.1 ) THEN >*/
if (kase == 1) {
/* Solve generalized Sylvester equation. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("N", &ijb, &n2, &n1, &a[i__ + i__ * a_dim1], lda,
&a[a_offset], lda, &work[1], &n2, &b[i__ + i__ *
b_dim1], ldb, &b[b_offset], ldb, &work[n1 * n2 +
1], &n2, &dscale, &dif[2], &work[(n1 << 1) * n2 +
1], &i__1, &iwork[1], &ierr, (ftnlen)1);
/*< ELSE >*/
} else {
/* Solve the transposed variant. */
/*< >*/
i__1 = *lwork - (n1 << 1) * n2;
dtgsyl_("T", &ijb, &n2, &n1, &a[i__ + i__ * a_dim1], lda,
&a[a_offset], lda, &work[1], &n2, &b[i__ + i__ *
b_dim1], ldb, &b[b_offset], ldb, &work[n1 * n2 +
1], &n2, &dscale, &dif[2], &work[(n1 << 1) * n2 +
1], &i__1, &iwork[1], &ierr, (ftnlen)1);
/*< END IF >*/
}
/*< GO TO 50 >*/
goto L50;
/*< END IF >*/
}
/*< DIF( 2 ) = DSCALE / DIF( 2 ) >*/
dif[2] = dscale / dif[2];
/*< END IF >*/
}
/*< END IF >*/
}
/*< 60 CONTINUE >*/
L60:
/* Compute generalized eigenvalues of reordered pair (A, B) and */
/* normalize the generalized Schur form. */
/*< PAIR = .FALSE. >*/
pair = FALSE_;
/*< DO 80 K = 1, N >*/
i__1 = *n;
for (k = 1; k <= i__1; ++k) {
/*< IF( PAIR ) THEN >*/
if (pair) {
/*< PAIR = .FALSE. >*/
pair = FALSE_;
/*< ELSE >*/
} else {
/*< IF( K.LT.N ) THEN >*/
if (k < *n) {
/*< IF( A( K+1, K ).NE.ZERO ) THEN >*/
if (a[k + 1 + k * a_dim1] != 0.) {
/*< PAIR = .TRUE. >*/
pair = TRUE_;
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( PAIR ) THEN >*/
if (pair) {
/* Compute the eigenvalue(s) at position K. */
/*< WORK( 1 ) = A( K, K ) >*/
work[1] = a[k + k * a_dim1];
/*< WORK( 2 ) = A( K+1, K ) >*/
work[2] = a[k + 1 + k * a_dim1];
/*< WORK( 3 ) = A( K, K+1 ) >*/
work[3] = a[k + (k + 1) * a_dim1];
/*< WORK( 4 ) = A( K+1, K+1 ) >*/
work[4] = a[k + 1 + (k + 1) * a_dim1];
/*< WORK( 5 ) = B( K, K ) >*/
work[5] = b[k + k * b_dim1];
/*< WORK( 6 ) = B( K+1, K ) >*/
work[6] = b[k + 1 + k * b_dim1];
/*< WORK( 7 ) = B( K, K+1 ) >*/
work[7] = b[k + (k + 1) * b_dim1];
/*< WORK( 8 ) = B( K+1, K+1 ) >*/
work[8] = b[k + 1 + (k + 1) * b_dim1];
/*< >*/
d__1 = smlnum * eps;
dlag2_(&work[1], &c__2, &work[5], &c__2, &d__1, &beta[k], &
beta[k + 1], &alphar[k], &alphar[k + 1], &alphai[k]);
/*< ALPHAI( K+1 ) = -ALPHAI( K ) >*/
alphai[k + 1] = -alphai[k];
/*< ELSE >*/
} else {
/*< IF( SIGN( ONE, B( K, K ) ).LT.ZERO ) THEN >*/
if (d_sign(&c_b28, &b[k + k * b_dim1]) < 0.) {
/* If B(K,K) is negative, make it positive */
/*< DO 70 I = 1, N >*/
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
/*< A( K, I ) = -A( K, I ) >*/
a[k + i__ * a_dim1] = -a[k + i__ * a_dim1];
/*< B( K, I ) = -B( K, I ) >*/
b[k + i__ * b_dim1] = -b[k + i__ * b_dim1];
/*< Q( I, K ) = -Q( I, K ) >*/
q[i__ + k * q_dim1] = -q[i__ + k * q_dim1];
/*< 70 CONTINUE >*/
/* L70: */
}
/*< END IF >*/
}
/*< ALPHAR( K ) = A( K, K ) >*/
alphar[k] = a[k + k * a_dim1];
/*< ALPHAI( K ) = ZERO >*/
alphai[k] = 0.;
/*< BETA( K ) = B( K, K ) >*/
beta[k] = b[k + k * b_dim1];
/*< END IF >*/
}
/*< END IF >*/
}
/*< 80 CONTINUE >*/
/* L80: */
}
/*< WORK( 1 ) = LWMIN >*/
work[1] = (doublereal) lwmin;
/*< IWORK( 1 ) = LIWMIN >*/
iwork[1] = liwmin;
/*< RETURN >*/
return 0;
/* End of DTGSEN */
/*< END >*/
} /* dtgsen_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -