📄 dtgsy2.c
字号:
/*< 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( I.GT.1 ) THEN >*/
if (i__ > 1) {
/*< >*/
i__2 = is - 1;
dger_(&i__2, &nb, &c_b27, &a[is * a_dim1 + 1], &c__1,
rhs, &c__1, &c__[js * c_dim1 + 1], ldc);
/*< >*/
i__2 = is - 1;
dger_(&i__2, &nb, &c_b27, &d__[is * d_dim1 + 1], &
c__1, rhs, &c__1, &f[js * f_dim1 + 1], ldf);
/*< END IF >*/
}
/*< IF( J.LT.Q ) THEN >*/
if (j < q) {
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[2], &b[js + (je + 1) * b_dim1],
ldb, &c__[is + (je + 1) * c_dim1], ldc);
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[2], &e[js + (je + 1) * e_dim1],
lde, &f[is + (je + 1) * f_dim1], ldf);
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[3], &b[jsp1 + (je + 1) * b_dim1],
ldb, &c__[is + (je + 1) * c_dim1], ldc);
/*< >*/
i__2 = *n - je;
daxpy_(&i__2, &rhs[3], &e[jsp1 + (je + 1) * e_dim1],
lde, &f[is + (je + 1) * f_dim1], ldf);
/*< 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( ISP1, IS ) >*/
z__[1] = a[isp1 + is * a_dim1];
/*< Z( 3, 1 ) = D( IS, IS ) >*/
z__[2] = d__[is + is * d_dim1];
/*< Z( 4, 1 ) = ZERO >*/
z__[3] = 0.;
/*< Z( 1, 2 ) = A( IS, ISP1 ) >*/
z__[8] = a[is + isp1 * a_dim1];
/*< Z( 2, 2 ) = A( ISP1, ISP1 ) >*/
z__[9] = a[isp1 + isp1 * a_dim1];
/*< Z( 3, 2 ) = D( IS, ISP1 ) >*/
z__[10] = d__[is + isp1 * d_dim1];
/*< Z( 4, 2 ) = D( ISP1, ISP1 ) >*/
z__[11] = d__[isp1 + isp1 * d_dim1];
/*< Z( 1, 3 ) = -B( JS, JS ) >*/
z__[16] = -b[js + js * b_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 ) = -B( JS, JS ) >*/
z__[25] = -b[js + js * b_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;
}
/*< IF( IJOB.EQ.0 ) THEN >*/
if (*ijob == 0) {
/*< >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 70 K = 1, N >*/
i__2 = *n;
for (k = 1; k <= i__2; ++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);
/*< 70 CONTINUE >*/
/* L70: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< >*/
dlatdf_(ijob, &zdim, z__, &c__8, rhs, rdsum, rdscal,
ipiv, jpiv);
/*< 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( I.GT.1 ) THEN >*/
if (i__ > 1) {
/*< >*/
i__2 = is - 1;
dgemv_("N", &i__2, &mb, &c_b27, &a[is * a_dim1 + 1],
lda, rhs, &c__1, &c_b42, &c__[js * c_dim1 + 1]
, &c__1, (ftnlen)1);
/*< >*/
i__2 = is - 1;
dgemv_("N", &i__2, &mb, &c_b27, &d__[is * d_dim1 + 1],
ldd, rhs, &c__1, &c_b42, &f[js * f_dim1 + 1],
&c__1, (ftnlen)1);
/*< END IF >*/
}
/*< IF( J.LT.Q ) THEN >*/
if (j < q) {
/*< >*/
i__2 = *n - je;
dger_(&mb, &i__2, &c_b42, &rhs[2], &c__1, &b[js + (je
+ 1) * b_dim1], ldb, &c__[is + (je + 1) *
c_dim1], ldc);
/*< >*/
i__2 = *n - je;
dger_(&mb, &i__2, &c_b42, &rhs[2], &c__1, &e[js + (je
+ 1) * e_dim1], ldb, &f[is + (je + 1) *
f_dim1], ldc);
/*< 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( ISP1, IS ) >*/
z__[1] = a[isp1 + is * a_dim1];
/*< Z( 5, 1 ) = D( IS, IS ) >*/
z__[4] = d__[is + is * d_dim1];
/*< Z( 1, 2 ) = A( IS, ISP1 ) >*/
z__[8] = a[is + isp1 * a_dim1];
/*< Z( 2, 2 ) = A( ISP1, ISP1 ) >*/
z__[9] = a[isp1 + isp1 * a_dim1];
/*< Z( 5, 2 ) = D( IS, ISP1 ) >*/
z__[12] = d__[is + isp1 * d_dim1];
/*< Z( 6, 2 ) = D( ISP1, ISP1 ) >*/
z__[13] = d__[isp1 + isp1 * d_dim1];
/*< Z( 3, 3 ) = A( IS, IS ) >*/
z__[18] = a[is + is * a_dim1];
/*< Z( 4, 3 ) = A( ISP1, IS ) >*/
z__[19] = a[isp1 + is * a_dim1];
/*< Z( 7, 3 ) = D( IS, IS ) >*/
z__[22] = d__[is + is * d_dim1];
/*< Z( 3, 4 ) = A( IS, ISP1 ) >*/
z__[26] = a[is + isp1 * a_dim1];
/*< Z( 4, 4 ) = A( ISP1, ISP1 ) >*/
z__[27] = a[isp1 + isp1 * a_dim1];
/*< Z( 7, 4 ) = D( IS, ISP1 ) >*/
z__[30] = d__[is + isp1 * d_dim1];
/*< Z( 8, 4 ) = D( ISP1, ISP1 ) >*/
z__[31] = d__[isp1 + isp1 * d_dim1];
/*< Z( 1, 5 ) = -B( JS, JS ) >*/
z__[32] = -b[js + js * b_dim1];
/*< Z( 3, 5 ) = -B( JS, JSP1 ) >*/
z__[34] = -b[js + jsp1 * b_dim1];
/*< Z( 5, 5 ) = -E( JS, JS ) >*/
z__[36] = -e[js + js * e_dim1];
/*< Z( 7, 5 ) = -E( JS, JSP1 ) >*/
z__[38] = -e[js + jsp1 * e_dim1];
/*< Z( 2, 6 ) = -B( JS, JS ) >*/
z__[41] = -b[js + js * b_dim1];
/*< Z( 4, 6 ) = -B( JS, JSP1 ) >*/
z__[43] = -b[js + jsp1 * b_dim1];
/*< Z( 6, 6 ) = -E( JS, JS ) >*/
z__[45] = -e[js + js * e_dim1];
/*< Z( 8, 6 ) = -E( JS, JSP1 ) >*/
z__[47] = -e[js + jsp1 * e_dim1];
/*< Z( 1, 7 ) = -B( JSP1, JS ) >*/
z__[48] = -b[jsp1 + js * b_dim1];
/*< Z( 3, 7 ) = -B( JSP1, JSP1 ) >*/
z__[50] = -b[jsp1 + jsp1 * b_dim1];
/*< Z( 7, 7 ) = -E( JSP1, JSP1 ) >*/
z__[54] = -e[jsp1 + jsp1 * e_dim1];
/*< Z( 2, 8 ) = -B( JSP1, JS ) >*/
z__[57] = -b[jsp1 + js * b_dim1];
/*< Z( 4, 8 ) = -B( JSP1, JSP1 ) >*/
z__[59] = -b[jsp1 + jsp1 * b_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 80 JJ = 0, NB - 1 >*/
i__2 = nb - 1;
for (jj = 0; jj <= i__2; ++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;
/*< 80 CONTINUE >*/
/* L80: */
}
/* 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;
}
/*< IF( IJOB.EQ.0 ) THEN >*/
if (*ijob == 0) {
/*< >*/
dgesc2_(&zdim, z__, &c__8, rhs, ipiv, jpiv, &scaloc);
/*< IF( SCALOC.NE.ONE ) THEN >*/
if (scaloc != 1.) {
/*< DO 90 K = 1, N >*/
i__2 = *n;
for (k = 1; k <= i__2; ++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);
/*< 90 CONTINUE >*/
/* L90: */
}
/*< SCALE = SCALE*SCALOC >*/
*scale *= scaloc;
/*< END IF >*/
}
/*< ELSE >*/
} else {
/*< >*/
dlatdf_(ijob, &zdim, z__, &c__8, rhs, rdsum, rdscal,
ipiv, jpiv);
/*< END IF >*/
}
/* Unpack solution vector(s) */
/*< K = 1 >*/
k = 1;
/*< II = MB*NB + 1 >*/
ii = mb * nb + 1;
/*< DO 100 JJ = 0, NB - 1 >*/
i__2 = nb - 1;
for (jj = 0; jj <= i__2; ++jj) {
/*< CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) >*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -