📄 dtgsy2.c
字号:
scaloc = 1.;
i__1 = p;
for (i = 1; i <= i__1; ++i) {
is = iwork[i];
isp1 = is + 1;
ie = iwork[i + 1] - 1;
mb = ie - is + 1;
i__2 = p + 2;
for (j = q; j >= i__2; --j) {
js = iwork[j];
jsp1 = js + 1;
je = iwork[j + 1] - 1;
nb = je - js + 1;
zdim = mb * nb << 1;
if (mb == 1 && nb == 1) {
/* Build a 2-by-2 system Z' * x = RHS */
z[0] = a[is + is * a_dim1];
z[1] = -b[js + js * b_dim1];
z[8] = d[is + is * d_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;
}
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__3 = *n;
for (k = 1; k <= i__3; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
/* 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 (j > p + 2) {
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[1];
i__3 = js - 1;
daxpy_(&i__3, &alpha, &e[js * e_dim1 + 1], &c__1, &f[is + f_dim1], ldf);
}
if (i < p) {
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[1];
i__3 = *m - ie;
daxpy_(&i__3, &alpha, &d[is + (ie + 1) * d_dim1], ldd, &c[ie + 1 + js * c_dim1], &c__1);
}
} 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] = -b[js + js * b_dim1];
z[3] = -b[jsp1 + js * b_dim1];
z[8] = 0.;
z[9] = a[is + is * a_dim1];
z[10] = -b[js + jsp1 * b_dim1];
z[11] = -b[jsp1 + jsp1 * b_dim1];
z[16] = d[is + is * d_dim1];
z[17] = 0.;
z[18] = -e[js + js * e_dim1];
z[19] = 0.;
z[24] = 0.;
z[25] = d[is + is * d_dim1];
z[26] = -e[js + jsp1 * e_dim1];
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;
}
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__3 = *n;
for (k = 1; k <= i__3; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
/* 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 (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);
}
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);
}
} 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[is + isp1 * a_dim1];
z[2] = -b[js + js * b_dim1];
z[3] = 0.;
z[8] = a[isp1 + is * a_dim1];
z[9] = a[isp1 + isp1 * a_dim1];
z[10] = 0.;
z[11] = -b[js + js * b_dim1];
z[16] = d[is + is * d_dim1];
z[17] = d[is + isp1 * d_dim1];
z[18] = -e[js + js * e_dim1];
z[19] = 0.;
z[24] = 0.;
z[25] = d[isp1 + isp1 * d_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;
}
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__3 = *n;
for (k = 1; k <= i__3; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
/* 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 (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);
}
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);
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);
}
} 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[is + isp1 * a_dim1];
z[4] = -b[js + js * b_dim1];
z[6] = -b[jsp1 + js * b_dim1];
z[8] = a[isp1 + is * a_dim1];
z[9] = a[isp1 + isp1 * a_dim1];
z[13] = -b[js + js * b_dim1];
z[15] = -b[jsp1 + js * b_dim1];
z[18] = a[is + is * a_dim1];
z[19] = a[is + isp1 * a_dim1];
z[20] = -b[js + jsp1 * b_dim1];
z[22] = -b[jsp1 + jsp1 * b_dim1];
z[26] = a[isp1 + is * a_dim1];
z[27] = a[isp1 + isp1 * a_dim1];
z[29] = -b[js + jsp1 * b_dim1];
z[31] = -b[jsp1 + jsp1 * b_dim1];
z[32] = d[is + is * d_dim1];
z[33] = d[is + isp1 * d_dim1];
z[36] = -e[js + js * e_dim1];
z[41] = d[isp1 + isp1 * d_dim1];
z[45] = -e[js + js * e_dim1];
z[50] = d[is + is * d_dim1];
z[51] = d[is + isp1 * d_dim1];
z[52] = -e[js + jsp1 * e_dim1];
z[54] = -e[jsp1 + jsp1 * e_dim1];
z[59] = d[isp1 + isp1 * d_dim1];
z[61] = -e[js + jsp1 * e_dim1];
z[63] = -e[jsp1 + jsp1 * e_dim1];
/* Set up right hand side(s) */
k = 1;
ii = mb * nb + 1;
i__3 = nb - 1;
for (jj = 0; jj <= i__3; ++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;
}
dgesc2_(&zdim, z, &c__8, rhs, ipiv, jpiv, &scaloc);
if (scaloc != 1.) {
i__3 = *n;
for (k = 1; k <= i__3; ++k) {
dscal_(m, &scaloc, &c[k * c_dim1 + 1], &c__1);
dscal_(m, &scaloc, &f[k * f_dim1 + 1], &c__1);
}
*scale *= scaloc;
}
/* Unpack solution vector(s) */
k = 1;
ii = mb * nb + 1;
i__3 = nb - 1;
for (jj = 0; jj <= i__3; ++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 (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);
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);
}
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);
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);
}
}
}
}
}
} /* dtgsy2_ */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -