📄 dtgsy2.c
字号:
dcopy_(&mb, &rhs[k - 1], &c__1, &c__[is + (js + jj) *
c_dim1], &c__1);
/*< CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) >*/
dcopy_(&mb, &rhs[ii - 1], &c__1, &f[is + (js + jj) *
f_dim1], &c__1);
/*< K = K + MB >*/
k += mb;
/*< II = II + MB >*/
ii += mb;
/*< 100 CONTINUE >*/
/* L100: */
}
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
/*< IF( I.GT.1 ) THEN >*/
if (i__ > 1) {
/*< >*/
i__2 = is - 1;
dgemm_("N", "N", &i__2, &nb, &mb, &c_b27, &a[is *
a_dim1 + 1], lda, rhs, &mb, &c_b42, &c__[js *
c_dim1 + 1], ldc, (ftnlen)1, (ftnlen)1);
/*< >*/
i__2 = is - 1;
dgemm_("N", "N", &i__2, &nb, &mb, &c_b27, &d__[is *
d_dim1 + 1], ldd, rhs, &mb, &c_b42, &f[js *
f_dim1 + 1], ldf, (ftnlen)1, (ftnlen)1);
/*< END IF >*/
}
/*< IF( J.LT.Q ) THEN >*/
if (j < q) {
/*< K = MB*NB + 1 >*/
k = mb * nb + 1;
/*< >*/
i__2 = *n - je;
dgemm_("N", "N", &mb, &i__2, &nb, &c_b42, &rhs[k - 1],
&mb, &b[js + (je + 1) * b_dim1], ldb, &c_b42,
&c__[is + (je + 1) * c_dim1], ldc, (ftnlen)1,
(ftnlen)1);
/*< >*/
i__2 = *n - je;
dgemm_("N", "N", &mb, &i__2, &nb, &c_b42, &rhs[k - 1],
&mb, &e[js + (je + 1) * e_dim1], lde, &c_b42,
&f[is + (je + 1) * f_dim1], ldf, (ftnlen)1, (
ftnlen)1);
/*< END IF >*/
}
/*< END IF >*/
}
/*< 110 CONTINUE >*/
/* L110: */
}
/*< 120 CONTINUE >*/
/* L120: */
}
/*< ELSE >*/
} else {
/* Solve (I, J) - subsystem */
/* A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J) */
/* R(I, I) * 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.;
/*< SCALOC = ONE >*/
scaloc = 1.;
/*< DO 200 I = 1, P >*/
i__1 = p;
for (i__ = 1; i__ <= i__1; ++i__) {
/*< IS = IWORK( I ) >*/
is = iwork[i__];
/*< ISP1 = IS + 1 >*/
isp1 = is + 1;
/*< IE = IWORK( I+1 ) - 1 >*/
ie = iwork[i__ + 1] - 1;
/*< MB = IE - IS + 1 >*/
mb = ie - is + 1;
/*< DO 190 J = Q, P + 2, -1 >*/
i__2 = p + 2;
for (j = q; j >= i__2; --j) {
/*< JS = IWORK( J ) >*/
js = iwork[j];
/*< JSP1 = JS + 1 >*/
jsp1 = js + 1;
/*< JE = IWORK( J+1 ) - 1 >*/
je = iwork[j + 1] - 1;
/*< NB = JE - JS + 1 >*/
nb = je - js + 1;
/*< ZDIM = MB*NB*2 >*/
zdim = mb * nb << 1;
/*< IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN >*/
if (mb == 1 && nb == 1) {
/* Build a 2-by-2 system Z' * x = RHS */
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = -B( JS, JS ) >*/
z__[1] = -b[js + js * b_dim1];
/*< Z( 1, 2 ) = D( IS, IS ) >*/
z__[8] = d__[is + is * d_dim1];
/*< Z( 2, 2 ) = -E( JS, JS ) >*/
z__[9] = -e[js + js * e_dim1];
/* Set up right hand side(s) */
/*< RHS( 1 ) = C( IS, JS ) >*/
rhs[0] = c__[is + js * c_dim1];
/*< RHS( 2 ) = F( IS, JS ) >*/
rhs[1] = f[is + js * f_dim1];
/* Solve Z' * x = RHS */
/*< CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) >*/
dgetc2_(&zdim, z__, &c__8, ipiv, jpiv, &ierr);
/*< >*/
if (ierr > 0) {
*info = ierr;
}
/*< CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 130 K = 1, N >*/
i__3 = *n;
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);
/*< 130 CONTINUE >*/
/* L130: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< C( IS, JS ) = RHS( 1 ) >*/
c__[is + js * c_dim1] = rhs[0];
/*< F( IS, JS ) = RHS( 2 ) >*/
f[is + js * f_dim1] = rhs[1];
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
/*< IF( J.GT.P+2 ) THEN >*/
if (j > p + 2) {
/*< ALPHA = RHS( 1 ) >*/
alpha = rhs[0];
/*< >*/
i__3 = js - 1;
daxpy_(&i__3, &alpha, &b[js * b_dim1 + 1], &c__1, &f[
is + f_dim1], ldf);
/*< ALPHA = RHS( 2 ) >*/
alpha = rhs[1];
/*< >*/
i__3 = js - 1;
daxpy_(&i__3, &alpha, &e[js * e_dim1 + 1], &c__1, &f[
is + f_dim1], ldf);
/*< END IF >*/
}
/*< IF( I.LT.P ) THEN >*/
if (i__ < p) {
/*< ALPHA = -RHS( 1 ) >*/
alpha = -rhs[0];
/*< >*/
i__3 = *m - ie;
daxpy_(&i__3, &alpha, &a[is + (ie + 1) * a_dim1], lda,
&c__[ie + 1 + js * c_dim1], &c__1);
/*< ALPHA = -RHS( 2 ) >*/
alpha = -rhs[1];
/*< >*/
i__3 = *m - ie;
daxpy_(&i__3, &alpha, &d__[is + (ie + 1) * d_dim1],
ldd, &c__[ie + 1 + js * c_dim1], &c__1);
/*< END IF >*/
}
/*< ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN >*/
} else if (mb == 1 && nb == 2) {
/* Build a 4-by-4 system Z' * x = RHS */
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = ZERO >*/
z__[1] = 0.;
/*< Z( 3, 1 ) = -B( JS, JS ) >*/
z__[2] = -b[js + js * b_dim1];
/*< Z( 4, 1 ) = -B( JSP1, JS ) >*/
z__[3] = -b[jsp1 + js * b_dim1];
/*< Z( 1, 2 ) = ZERO >*/
z__[8] = 0.;
/*< Z( 2, 2 ) = A( IS, IS ) >*/
z__[9] = a[is + is * a_dim1];
/*< Z( 3, 2 ) = -B( JS, JSP1 ) >*/
z__[10] = -b[js + jsp1 * b_dim1];
/*< Z( 4, 2 ) = -B( JSP1, JSP1 ) >*/
z__[11] = -b[jsp1 + jsp1 * b_dim1];
/*< Z( 1, 3 ) = D( IS, IS ) >*/
z__[16] = d__[is + is * d_dim1];
/*< Z( 2, 3 ) = ZERO >*/
z__[17] = 0.;
/*< Z( 3, 3 ) = -E( JS, JS ) >*/
z__[18] = -e[js + js * e_dim1];
/*< Z( 4, 3 ) = ZERO >*/
z__[19] = 0.;
/*< Z( 1, 4 ) = ZERO >*/
z__[24] = 0.;
/*< Z( 2, 4 ) = D( IS, IS ) >*/
z__[25] = d__[is + is * d_dim1];
/*< Z( 3, 4 ) = -E( JS, JSP1 ) >*/
z__[26] = -e[js + jsp1 * e_dim1];
/*< Z( 4, 4 ) = -E( JSP1, JSP1 ) >*/
z__[27] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
/*< RHS( 1 ) = C( IS, JS ) >*/
rhs[0] = c__[is + js * c_dim1];
/*< RHS( 2 ) = C( IS, JSP1 ) >*/
rhs[1] = c__[is + jsp1 * c_dim1];
/*< RHS( 3 ) = F( IS, JS ) >*/
rhs[2] = f[is + js * f_dim1];
/*< RHS( 4 ) = F( IS, JSP1 ) >*/
rhs[3] = f[is + jsp1 * f_dim1];
/* Solve Z' * x = RHS */
/*< CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) >*/
dgetc2_(&zdim, z__, &c__8, ipiv, jpiv, &ierr);
/*< >*/
if (ierr > 0) {
*info = ierr;
}
/*< CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 140 K = 1, N >*/
i__3 = *n;
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);
/*< 140 CONTINUE >*/
/* L140: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< C( IS, JS ) = RHS( 1 ) >*/
c__[is + js * c_dim1] = rhs[0];
/*< C( IS, JSP1 ) = RHS( 2 ) >*/
c__[is + jsp1 * c_dim1] = rhs[1];
/*< F( IS, JS ) = RHS( 3 ) >*/
f[is + js * f_dim1] = rhs[2];
/*< F( IS, JSP1 ) = RHS( 4 ) >*/
f[is + jsp1 * f_dim1] = rhs[3];
/* 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;
daxpy_(&i__3, rhs, &b[js * b_dim1 + 1], &c__1, &f[is
+ f_dim1], ldf);
/*< >*/
i__3 = js - 1;
daxpy_(&i__3, &rhs[1], &b[jsp1 * b_dim1 + 1], &c__1, &
f[is + f_dim1], ldf);
/*< >*/
i__3 = js - 1;
daxpy_(&i__3, &rhs[2], &e[js * e_dim1 + 1], &c__1, &f[
is + f_dim1], ldf);
/*< >*/
i__3 = js - 1;
daxpy_(&i__3, &rhs[3], &e[jsp1 * e_dim1 + 1], &c__1, &
f[is + f_dim1], ldf);
/*< END IF >*/
}
/*< IF( I.LT.P ) THEN >*/
if (i__ < p) {
/*< >*/
i__3 = *m - ie;
dger_(&i__3, &nb, &c_b27, &a[is + (ie + 1) * a_dim1],
lda, rhs, &c__1, &c__[ie + 1 + js * c_dim1],
ldc);
/*< >*/
i__3 = *m - ie;
dger_(&i__3, &nb, &c_b27, &d__[is + (ie + 1) * d_dim1]
, ldd, &rhs[2], &c__1, &c__[ie + 1 + js *
c_dim1], ldc);
/*< END IF >*/
}
/*< ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN >*/
} else if (mb == 2 && nb == 1) {
/* Build a 4-by-4 system Z' * x = RHS */
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = A( IS, ISP1 ) >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -