📄 dtgsyl.c
字号:
/*< PPQQ = 0 >*/
ppqq = 0;
/*< >*/
dtgsy2_(trans, &ifunc, &mb, &nb, &a[is + is * a_dim1],
lda, &b[js + js * b_dim1], ldb, &c__[is + js *
c_dim1], ldc, &d__[is + is * d_dim1], ldd, &e[js
+ js * e_dim1], lde, &f[is + js * f_dim1], ldf, &
scaloc, &dsum, &dscale, &iwork[q + 2], &ppqq, &
linfo, (ftnlen)1);
/*< >*/
if (linfo > 0) {
*info = linfo;
}
/*< PQ = PQ + PPQQ >*/
pq += ppqq;
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 80 K = 1, JS - 1 >*/
i__3 = js - 1;
for (k = 1; k <= i__3; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 80 CONTINUE >*/
/* L80: */
}
/*< DO 90 K = JS, JE >*/
i__3 = je;
for (k = js; k <= i__3; ++k) {
/*< CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 ) >*/
i__4 = is - 1;
dscal_(&i__4, &scaloc, &c__[k * c_dim1 + 1], &
c__1);
/*< CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 ) >*/
i__4 = is - 1;
dscal_(&i__4, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 90 CONTINUE >*/
/* L90: */
}
/*< DO 100 K = JS, JE >*/
i__3 = je;
for (k = js; k <= i__3; ++k) {
/*< CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 ) >*/
i__4 = *m - ie;
dscal_(&i__4, &scaloc, &c__[ie + 1 + k * c_dim1],
&c__1);
/*< CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 ) >*/
i__4 = *m - ie;
dscal_(&i__4, &scaloc, &f[ie + 1 + k * f_dim1], &
c__1);
/*< 100 CONTINUE >*/
/* L100: */
}
/*< DO 110 K = JE + 1, N >*/
i__3 = *n;
for (k = je + 1; k <= i__3; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 110 CONTINUE >*/
/* L110: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
/*< IF( I.GT.1 ) THEN >*/
if (i__ > 1) {
/*< >*/
i__3 = is - 1;
dgemm_("N", "N", &i__3, &nb, &mb, &c_b53, &a[is *
a_dim1 + 1], lda, &c__[is + js * c_dim1], ldc,
&c_b54, &c__[js * c_dim1 + 1], ldc, (ftnlen)
1, (ftnlen)1);
/*< >*/
i__3 = is - 1;
dgemm_("N", "N", &i__3, &nb, &mb, &c_b53, &d__[is *
d_dim1 + 1], ldd, &c__[is + js * c_dim1], ldc,
&c_b54, &f[js * f_dim1 + 1], ldf, (ftnlen)1,
(ftnlen)1);
/*< END IF >*/
}
/*< IF( J.LT.Q ) THEN >*/
if (j < q) {
/*< >*/
i__3 = *n - je;
dgemm_("N", "N", &mb, &i__3, &nb, &c_b54, &f[is + js *
f_dim1], ldf, &b[js + (je + 1) * b_dim1],
ldb, &c_b54, &c__[is + (je + 1) * c_dim1],
ldc, (ftnlen)1, (ftnlen)1);
/*< >*/
i__3 = *n - je;
dgemm_("N", "N", &mb, &i__3, &nb, &c_b54, &f[is + js *
f_dim1], ldf, &e[js + (je + 1) * e_dim1],
lde, &c_b54, &f[is + (je + 1) * f_dim1], ldf,
(ftnlen)1, (ftnlen)1);
/*< END IF >*/
}
/*< 120 CONTINUE >*/
/* L120: */
}
/*< 130 CONTINUE >*/
/* L130: */
}
/*< IF( DSCALE.NE.ZERO ) THEN >*/
if (dscale != 0.) {
/*< IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN >*/
if (*ijob == 1 || *ijob == 3) {
/*< DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) >*/
*dif = sqrt((doublereal) ((*m << 1) * *n)) / (dscale *
sqrt(dsum));
/*< ELSE >*/
} else {
/*< DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) >*/
*dif = sqrt((doublereal) pq) / (dscale * sqrt(dsum));
/*< END IF >*/
}
/*< END IF >*/
}
/*< IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN >*/
if (isolve == 2 && iround == 1) {
/*< IFUNC = IJOB >*/
ifunc = *ijob;
/*< SCALE2 = SCALE >*/
scale2 = *scale;
/*< CALL DLACPY( 'F', M, N, C, LDC, WORK, M ) >*/
dlacpy_("F", m, n, &c__[c_offset], ldc, &work[1], m, (ftnlen)
1);
/*< CALL DLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) >*/
dlacpy_("F", m, n, &f[f_offset], ldf, &work[*m * *n + 1], m, (
ftnlen)1);
/*< DO 140 J = 1, N >*/
i__2 = *n;
for (j = 1; j <= i__2; ++j) {
/*< CALL DCOPY( M, ZERO, 0, C( 1, J ), 1 ) >*/
dcopy_(m, &c_b14, &c__0, &c__[j * c_dim1 + 1], &c__1);
/*< CALL DCOPY( M, ZERO, 0, F( 1, J ), 1 ) >*/
dcopy_(m, &c_b14, &c__0, &f[j * f_dim1 + 1], &c__1);
/*< 140 CONTINUE >*/
/* L140: */
}
/*< ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN >*/
} else if (isolve == 2 && iround == 2) {
/*< CALL DLACPY( 'F', M, N, WORK, M, C, LDC ) >*/
dlacpy_("F", m, n, &work[1], m, &c__[c_offset], ldc, (ftnlen)
1);
/*< CALL DLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) >*/
dlacpy_("F", m, n, &work[*m * *n + 1], m, &f[f_offset], ldf, (
ftnlen)1);
/*< SCALE = SCALE2 >*/
*scale = scale2;
/*< END IF >*/
}
/*< 150 CONTINUE >*/
/* L150: */
}
/*< ELSE >*/
} else {
/* Solve transposed (I, J)-subsystem */
/* A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J) */
/* R(I, J) * B(J, J)' + L(I, J) * E(J, J)' = -F(I, J) */
/* for I = 1,2,..., P; J = Q, Q-1,..., 1 */
/*< SCALE = ONE >*/
*scale = 1.;
/*< DO 210 I = 1, P >*/
i__1 = p;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< IS = IWORK( I ) >*/
is = iwork[i__];
/*< IE = IWORK( I+1 ) - 1 >*/
ie = iwork[i__ + 1] - 1;
/*< MB = IE - IS + 1 >*/
mb = ie - is + 1;
/*< DO 200 J = Q, P + 2, -1 >*/
i__2 = p + 2;
for (j = q; j >= i__2; --j) {
/*< JS = IWORK( J ) >*/
js = iwork[j];
/*< JE = IWORK( J+1 ) - 1 >*/
je = iwork[j + 1] - 1;
/*< NB = JE - JS + 1 >*/
nb = je - js + 1;
/*< >*/
dtgsy2_(trans, &ifunc, &mb, &nb, &a[is + is * a_dim1], lda, &
b[js + js * b_dim1], ldb, &c__[is + js * c_dim1], ldc,
&d__[is + is * d_dim1], ldd, &e[js + js * e_dim1],
lde, &f[is + js * f_dim1], ldf, &scaloc, &dsum, &
dscale, &iwork[q + 2], &ppqq, &linfo, (ftnlen)1);
/*< >*/
if (linfo > 0) {
*info = linfo;
}
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 160 K = 1, JS - 1 >*/
i__3 = js - 1;
for (k = 1; k <= i__3; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 160 CONTINUE >*/
/* L160: */
}
/*< DO 170 K = JS, JE >*/
i__3 = je;
for (k = js; k <= i__3; ++k) {
/*< CALL DSCAL( IS-1, SCALOC, C( 1, K ), 1 ) >*/
i__4 = is - 1;
dscal_(&i__4, &scaloc, &c__[k * c_dim1 + 1], &c__1);
/*< CALL DSCAL( IS-1, SCALOC, F( 1, K ), 1 ) >*/
i__4 = is - 1;
dscal_(&i__4, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 170 CONTINUE >*/
/* L170: */
}
/*< DO 180 K = JS, JE >*/
i__3 = je;
for (k = js; k <= i__3; ++k) {
/*< CALL DSCAL( M-IE, SCALOC, C( IE+1, K ), 1 ) >*/
i__4 = *m - ie;
dscal_(&i__4, &scaloc, &c__[ie + 1 + k * c_dim1], &
c__1);
/*< CALL DSCAL( M-IE, SCALOC, F( IE+1, K ), 1 ) >*/
i__4 = *m - ie;
dscal_(&i__4, &scaloc, &f[ie + 1 + k * f_dim1], &c__1)
;
/*< 180 CONTINUE >*/
/* L180: */
}
/*< DO 190 K = JE + 1, N >*/
i__3 = *n;
for (k = je + 1; k <= i__3; ++k) {
/*< CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &c__[k * c_dim1 + 1], &c__1);
/*< CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) >*/
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
/*< 190 CONTINUE >*/
/* L190: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Substitute R(I, J) and L(I, J) into remaining equation. */
/*< IF( J.GT.P+2 ) THEN >*/
if (j > p + 2) {
/*< >*/
i__3 = js - 1;
dgemm_("N", "T", &mb, &i__3, &nb, &c_b54, &c__[is + js *
c_dim1], ldc, &b[js * b_dim1 + 1], ldb, &c_b54, &
f[is + f_dim1], ldf, (ftnlen)1, (ftnlen)1);
/*< >*/
i__3 = js - 1;
dgemm_("N", "T", &mb, &i__3, &nb, &c_b54, &f[is + js *
f_dim1], ldf, &e[js * e_dim1 + 1], lde, &c_b54, &
f[is + f_dim1], ldf, (ftnlen)1, (ftnlen)1);
/*< END IF >*/
}
/*< IF( I.LT.P ) THEN >*/
if (i__ < p) {
/*< >*/
i__3 = *m - ie;
dgemm_("T", "N", &i__3, &nb, &mb, &c_b53, &a[is + (ie + 1)
* a_dim1], lda, &c__[is + js * c_dim1], ldc, &
c_b54, &c__[ie + 1 + js * c_dim1], ldc, (ftnlen)1,
(ftnlen)1);
/*< >*/
i__3 = *m - ie;
dgemm_("T", "N", &i__3, &nb, &mb, &c_b53, &d__[is + (ie +
1) * d_dim1], ldd, &f[is + js * f_dim1], ldf, &
c_b54, &c__[ie + 1 + js * c_dim1], ldc, (ftnlen)1,
(ftnlen)1);
/*< END IF >*/
}
/*< 200 CONTINUE >*/
/* L200: */
}
/*< 210 CONTINUE >*/
/* L210: */
}
/*< END IF >*/
}
/*< WORK( 1 ) = LWMIN >*/
work[1] = (doublereal) lwmin;
/*< RETURN >*/
return 0;
/* End of DTGSYL */
/*< END >*/
} /* dtgsyl_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -