📄 dtgsy2.c
字号:
z__[1] = a[is + isp1 * a_dim1];
/*< Z( 3, 1 ) = -B( JS, JS ) >*/
z__[2] = -b[js + js * b_dim1];
/*< Z( 4, 1 ) = ZERO >*/
z__[3] = 0.;
/*< Z( 1, 2 ) = A( ISP1, IS ) >*/
z__[8] = a[isp1 + is * a_dim1];
/*< Z( 2, 2 ) = A( ISP1, ISP1 ) >*/
z__[9] = a[isp1 + isp1 * a_dim1];
/*< Z( 3, 2 ) = ZERO >*/
z__[10] = 0.;
/*< Z( 4, 2 ) = -B( JS, JS ) >*/
z__[11] = -b[js + js * b_dim1];
/*< Z( 1, 3 ) = D( IS, IS ) >*/
z__[16] = d__[is + is * d_dim1];
/*< Z( 2, 3 ) = D( IS, ISP1 ) >*/
z__[17] = d__[is + isp1 * d_dim1];
/*< 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( ISP1, ISP1 ) >*/
z__[25] = d__[isp1 + isp1 * d_dim1];
/*< Z( 3, 4 ) = ZERO >*/
z__[26] = 0.;
/*< Z( 4, 4 ) = -E( JS, JS ) >*/
z__[27] = -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 ) = C( ISP1, JS ) >*/
rhs[1] = c__[isp1 + js * c_dim1];
/*< RHS( 3 ) = F( IS, JS ) >*/
rhs[2] = f[is + js * f_dim1];
/*< RHS( 4 ) = F( ISP1, JS ) >*/
rhs[3] = f[isp1 + 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 150 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);
/*< 150 CONTINUE >*/
/* L150: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< C( IS, JS ) = RHS( 1 ) >*/
c__[is + js * c_dim1] = rhs[0];
/*< C( ISP1, JS ) = RHS( 2 ) >*/
c__[isp1 + js * c_dim1] = rhs[1];
/*< F( IS, JS ) = RHS( 3 ) >*/
f[is + js * f_dim1] = rhs[2];
/*< F( ISP1, JS ) = RHS( 4 ) >*/
f[isp1 + js * 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;
dger_(&mb, &i__3, &c_b42, rhs, &c__1, &b[js * b_dim1
+ 1], &c__1, &f[is + f_dim1], ldf);
/*< >*/
i__3 = js - 1;
dger_(&mb, &i__3, &c_b42, &rhs[2], &c__1, &e[js *
e_dim1 + 1], &c__1, &f[is + f_dim1], ldf);
/*< END IF >*/
}
/*< IF( I.LT.P ) THEN >*/
if (i__ < p) {
/*< >*/
i__3 = *m - ie;
dgemv_("T", &mb, &i__3, &c_b27, &a[is + (ie + 1) *
a_dim1], lda, rhs, &c__1, &c_b42, &c__[ie + 1
+ js * c_dim1], &c__1, (ftnlen)1);
/*< >*/
i__3 = *m - ie;
dgemv_("T", &mb, &i__3, &c_b27, &d__[is + (ie + 1) *
d_dim1], ldd, &rhs[2], &c__1, &c_b42, &c__[ie
+ 1 + js * c_dim1], &c__1, (ftnlen)1);
/*< END IF >*/
}
/*< ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN >*/
} else if (mb == 2 && nb == 2) {
/* Build an 8-by-8 system Z' * x = RHS */
/*< CALL DCOPY( LDZ*LDZ, ZERO, 0, Z, 1 ) >*/
dcopy_(&c__64, &c_b54, &c__0, z__, &c__1);
/*< Z( 1, 1 ) = A( IS, IS ) >*/
z__[0] = a[is + is * a_dim1];
/*< Z( 2, 1 ) = A( IS, ISP1 ) >*/
z__[1] = a[is + isp1 * a_dim1];
/*< Z( 5, 1 ) = -B( JS, JS ) >*/
z__[4] = -b[js + js * b_dim1];
/*< Z( 7, 1 ) = -B( JSP1, JS ) >*/
z__[6] = -b[jsp1 + js * b_dim1];
/*< Z( 1, 2 ) = A( ISP1, IS ) >*/
z__[8] = a[isp1 + is * a_dim1];
/*< Z( 2, 2 ) = A( ISP1, ISP1 ) >*/
z__[9] = a[isp1 + isp1 * a_dim1];
/*< Z( 6, 2 ) = -B( JS, JS ) >*/
z__[13] = -b[js + js * b_dim1];
/*< Z( 8, 2 ) = -B( JSP1, JS ) >*/
z__[15] = -b[jsp1 + js * b_dim1];
/*< Z( 3, 3 ) = A( IS, IS ) >*/
z__[18] = a[is + is * a_dim1];
/*< Z( 4, 3 ) = A( IS, ISP1 ) >*/
z__[19] = a[is + isp1 * a_dim1];
/*< Z( 5, 3 ) = -B( JS, JSP1 ) >*/
z__[20] = -b[js + jsp1 * b_dim1];
/*< Z( 7, 3 ) = -B( JSP1, JSP1 ) >*/
z__[22] = -b[jsp1 + jsp1 * b_dim1];
/*< Z( 3, 4 ) = A( ISP1, IS ) >*/
z__[26] = a[isp1 + is * a_dim1];
/*< Z( 4, 4 ) = A( ISP1, ISP1 ) >*/
z__[27] = a[isp1 + isp1 * a_dim1];
/*< Z( 6, 4 ) = -B( JS, JSP1 ) >*/
z__[29] = -b[js + jsp1 * b_dim1];
/*< Z( 8, 4 ) = -B( JSP1, JSP1 ) >*/
z__[31] = -b[jsp1 + jsp1 * b_dim1];
/*< Z( 1, 5 ) = D( IS, IS ) >*/
z__[32] = d__[is + is * d_dim1];
/*< Z( 2, 5 ) = D( IS, ISP1 ) >*/
z__[33] = d__[is + isp1 * d_dim1];
/*< Z( 5, 5 ) = -E( JS, JS ) >*/
z__[36] = -e[js + js * e_dim1];
/*< Z( 2, 6 ) = D( ISP1, ISP1 ) >*/
z__[41] = d__[isp1 + isp1 * d_dim1];
/*< Z( 6, 6 ) = -E( JS, JS ) >*/
z__[45] = -e[js + js * e_dim1];
/*< Z( 3, 7 ) = D( IS, IS ) >*/
z__[50] = d__[is + is * d_dim1];
/*< Z( 4, 7 ) = D( IS, ISP1 ) >*/
z__[51] = d__[is + isp1 * d_dim1];
/*< Z( 5, 7 ) = -E( JS, JSP1 ) >*/
z__[52] = -e[js + jsp1 * e_dim1];
/*< Z( 7, 7 ) = -E( JSP1, JSP1 ) >*/
z__[54] = -e[jsp1 + jsp1 * e_dim1];
/*< Z( 4, 8 ) = D( ISP1, ISP1 ) >*/
z__[59] = d__[isp1 + isp1 * d_dim1];
/*< Z( 6, 8 ) = -E( JS, JSP1 ) >*/
z__[61] = -e[js + jsp1 * e_dim1];
/*< Z( 8, 8 ) = -E( JSP1, JSP1 ) >*/
z__[63] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
/*< K = 1 >*/
k = 1;
/*< II = MB*NB + 1 >*/
ii = mb * nb + 1;
/*< DO 160 JJ = 0, NB - 1 >*/
i__3 = nb - 1;
for (jj = 0; jj <= i__3; ++jj) {
/*< CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) >*/
dcopy_(&mb, &c__[is + (js + jj) * c_dim1], &c__1, &
rhs[k - 1], &c__1);
/*< CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) >*/
dcopy_(&mb, &f[is + (js + jj) * f_dim1], &c__1, &rhs[
ii - 1], &c__1);
/*< K = K + MB >*/
k += mb;
/*< II = II + MB >*/
ii += mb;
/*< 160 CONTINUE >*/
/* L160: */
}
/* 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 170 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);
/*< 170 CONTINUE >*/
/* L170: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< K = 1 >*/
k = 1;
/*< II = MB*NB + 1 >*/
ii = mb * nb + 1;
/*< DO 180 JJ = 0, NB - 1 >*/
i__3 = nb - 1;
for (jj = 0; jj <= i__3; ++jj) {
/*< CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) >*/
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;
/*< 180 CONTINUE >*/
/* L180: */
}
/* 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_b42, &c__[is +
js * c_dim1], ldc, &b[js * b_dim1 + 1], ldb, &
c_b42, &f[is + f_dim1], ldf, (ftnlen)1, (
ftnlen)1);
/*< >*/
i__3 = js - 1;
dgemm_("N", "T", &mb, &i__3, &nb, &c_b42, &f[is + js *
f_dim1], ldf, &e[js * e_dim1 + 1], lde, &
c_b42, &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_b27, &a[is + (ie
+ 1) * a_dim1], lda, &c__[is + js * c_dim1],
ldc, &c_b42, &c__[ie + 1 + js * c_dim1], ldc,
(ftnlen)1, (ftnlen)1);
/*< >*/
i__3 = *m - ie;
dgemm_("T", "N", &i__3, &nb, &mb, &c_b27, &d__[is + (
ie + 1) * d_dim1], ldd, &f[is + js * f_dim1],
ldf, &c_b42, &c__[ie + 1 + js * c_dim1], ldc,
(ftnlen)1, (ftnlen)1);
/*< END IF >*/
}
/*< END IF >*/
}
/*< 190 CONTINUE >*/
/* L190: */
}
/*< 200 CONTINUE >*/
/* L200: */
}
/*< END IF >*/
}
/*< RETURN >*/
return 0;
/* End of DTGSY2 */
/*< END >*/
} /* dtgsy2_ */
#ifdef __cplusplus
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -