📄 dtgsy2.c
字号:
/* Build a 2-by-2 system Z * x = RHS */
z[0] = a[is + is * a_dim1];
z[1] = d[is + is * d_dim1];
z[8] = -b[js + js * b_dim1];
z[9] = -e[js + js * e_dim1];
/* Set up right hand side(s) */
rhs[0] = c[is + js * c_dim1];
rhs[1] = f[is + js * f_dim1];
/* Solve Z * x = RHS */
dgetc2_(&zdim, z, &c__8, ipiv, jpiv, &ierr);
if (ierr > 0) {
*info = ierr;
}
if (*ijob == 0) {
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
} else {
dlatdf_(ijob, &zdim, z, &c__8, rhs, rdsum, rdscal, ipiv, jpiv);
}
/* Unpack solution vector(s) */
c[is + js * c_dim1] = rhs[0];
f[is + js * f_dim1] = rhs[1];
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
if (i > 1) {
alpha = -rhs[0];
i__2 = is - 1;
daxpy_(&i__2, &alpha, &a[is * a_dim1 + 1], &c__1, &c[js * c_dim1 + 1], &c__1);
i__2 = is - 1;
daxpy_(&i__2, &alpha, &d[is * d_dim1 + 1], &c__1, &f[js * f_dim1 + 1], &c__1);
}
if (j < q) {
i__2 = *n - je;
daxpy_(&i__2, &rhs[1], &b[js + (je + 1) * b_dim1], ldb, &c[is + (je + 1) * c_dim1], ldc);
i__2 = *n - je;
daxpy_(&i__2, &rhs[1], &e[js + (je + 1) * e_dim1], lde, &f[is + (je + 1) * f_dim1], ldf);
}
} else if (mb == 1 && nb == 2) {
/* Build a 4-by-4 system Z * x = RHS */
z[0] = a[is + is * a_dim1];
z[1] = 0.;
z[2] = d[is + is * d_dim1];
z[3] = 0.;
z[8] = 0.;
z[9] = a[is + is * a_dim1];
z[10] = 0.;
z[11] = d[is + is * d_dim1];
z[16] = -b[js + js * b_dim1];
z[17] = -b[js + jsp1 * b_dim1];
z[18] = -e[js + js * e_dim1];
z[19] = -e[js + jsp1 * e_dim1];
z[24] = -b[jsp1 + js * b_dim1];
z[25] = -b[jsp1 + jsp1 * b_dim1];
z[26] = 0.;
z[27] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
rhs[0] = c[is + js * c_dim1];
rhs[1] = c[is + jsp1 * c_dim1];
rhs[2] = f[is + js * f_dim1];
rhs[3] = f[is + jsp1 * f_dim1];
/* Solve Z * x = RHS */
dgetc2_(&zdim, z, &c__8, ipiv, jpiv, &ierr);
if (ierr > 0) {
*info = ierr;
}
if (*ijob == 0) {
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
} else {
dlatdf_(ijob, &zdim, z, &c__8, rhs, rdsum, rdscal, ipiv, jpiv);
}
/* Unpack solution vector(s) */
c[is + js * c_dim1] = rhs[0];
c[is + jsp1 * c_dim1] = rhs[1];
f[is + js * f_dim1] = rhs[2];
f[is + jsp1 * f_dim1] = rhs[3];
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
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);
}
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);
}
} else if (mb == 2 && nb == 1) {
/* Build a 4-by-4 system Z * x = RHS */
z[0] = a[is + is * a_dim1];
z[1] = a[isp1 + is * a_dim1];
z[2] = d[is + is * d_dim1];
z[3] = 0.;
z[8] = a[is + isp1 * a_dim1];
z[9] = a[isp1 + isp1 * a_dim1];
z[10] = d[is + isp1 * d_dim1];
z[11] = d[isp1 + isp1 * d_dim1];
z[16] = -b[js + js * b_dim1];
z[17] = 0.;
z[18] = -e[js + js * e_dim1];
z[19] = 0.;
z[24] = 0.;
z[25] = -b[js + js * b_dim1];
z[26] = 0.;
z[27] = -e[js + js * e_dim1];
/* Set up right hand side(s) */
rhs[0] = c[is + js * c_dim1];
rhs[1] = c[isp1 + js * c_dim1];
rhs[2] = f[is + js * f_dim1];
rhs[3] = f[isp1 + js * f_dim1];
/* Solve Z * x = RHS */
dgetc2_(&zdim, z, &c__8, ipiv, jpiv, &ierr);
if (ierr > 0) {
*info = ierr;
}
if (*ijob == 0) {
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
} else {
dlatdf_(ijob, &zdim, z, &c__8, rhs, rdsum, rdscal, ipiv, jpiv);
}
/* Unpack solution vector(s) */
c[is + js * c_dim1] = rhs[0];
c[isp1 + js * c_dim1] = rhs[1];
f[is + js * f_dim1] = rhs[2];
f[isp1 + js * f_dim1] = rhs[3];
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
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);
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);
}
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);
}
} else if (mb == 2 && nb == 2) {
/* Build an 8-by-8 system Z * x = RHS */
dcopy_(&c__64, &c_b54, &c__0, z, &c__1);
z[0] = a[is + is * a_dim1];
z[1] = a[isp1 + is * a_dim1];
z[4] = d[is + is * d_dim1];
z[8] = a[is + isp1 * a_dim1];
z[9] = a[isp1 + isp1 * a_dim1];
z[12] = d[is + isp1 * d_dim1];
z[13] = d[isp1 + isp1 * d_dim1];
z[18] = a[is + is * a_dim1];
z[19] = a[isp1 + is * a_dim1];
z[22] = d[is + is * d_dim1];
z[26] = a[is + isp1 * a_dim1];
z[27] = a[isp1 + isp1 * a_dim1];
z[30] = d[is + isp1 * d_dim1];
z[31] = d[isp1 + isp1 * d_dim1];
z[32] = -b[js + js * b_dim1];
z[34] = -b[js + jsp1 * b_dim1];
z[36] = -e[js + js * e_dim1];
z[38] = -e[js + jsp1 * e_dim1];
z[41] = -b[js + js * b_dim1];
z[43] = -b[js + jsp1 * b_dim1];
z[45] = -e[js + js * e_dim1];
z[47] = -e[js + jsp1 * e_dim1];
z[48] = -b[jsp1 + js * b_dim1];
z[50] = -b[jsp1 + jsp1 * b_dim1];
z[54] = -e[jsp1 + jsp1 * e_dim1];
z[57] = -b[jsp1 + js * b_dim1];
z[59] = -b[jsp1 + jsp1 * b_dim1];
z[63] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
k = 1;
ii = mb * nb + 1;
i__2 = nb - 1;
for (jj = 0; jj <= i__2; ++jj) {
dcopy_(&mb, &c[is + (js + jj) * c_dim1], &c__1, &rhs[k - 1], &c__1);
dcopy_(&mb, &f[is + (js + jj) * f_dim1], &c__1, &rhs[ii - 1], &c__1);
k += mb;
ii += mb;
}
/* Solve Z * x = RHS */
dgetc2_(&zdim, z, &c__8, ipiv, jpiv, &ierr);
if (ierr > 0) {
*info = ierr;
}
if (*ijob == 0) {
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
} else {
dlatdf_(ijob, &zdim, z, &c__8, rhs, rdsum, rdscal, ipiv, jpiv);
}
/* Unpack solution vector(s) */
k = 1;
ii = mb * nb + 1;
i__2 = nb - 1;
for (jj = 0; jj <= i__2; ++jj) {
dcopy_(&mb, &rhs[k - 1], &c__1, &c[is + (js + jj) * c_dim1], &c__1);
dcopy_(&mb, &rhs[ii - 1], &c__1, &f[is + (js + jj) * f_dim1], &c__1);
k += mb;
ii += mb;
}
/* Substitute R(I, J) and L(I, J) into remaining */
/* equation. */
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);
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);
}
if (j < q) {
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);
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);
}
}
}
}
} 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 = 1.;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -