dlarfb.c
来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 553 行 · 第 1/2 页
C
553 行
/* C2 := C2 - W' */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *n; ++i) {
c[*m - *k + j + i * *ldc] -= work[i + j * *ldwork];
}
}
} else if (lsame_(side, "R")) {
/* Form C * H or C * H' where C = ( C1 C2 ) */
/* W := C * V = (C1*V1 + C2*V2) (stored in WORK) */
/* W := C2 */
for (j = 0; j < *k; ++j) {
dcopy_(m, &c[(*n - *k + j) * *ldc], &c__1, &work[j * *ldwork], &c__1);
}
/* W := W * V2 */
dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14, &v[*n - *k], ldv, work, ldwork);
if (*n > *k) {
/* W := W + C1 * V1 */
i__1 = *n - *k;
dgemm_("No transpose", "No transpose", m, k, &i__1, &c_b14, c, ldc, v, ldv, &c_b14, work, ldwork);
}
/* W := W * T or W * T' */
dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, t, ldt, work, ldwork);
/* C := C - W * V' */
if (*n > *k) {
/* C1 := C1 - W * V1' */
i__1 = *n - *k;
dgemm_("No transpose", "Transpose", m, &i__1, k, &c_b25, work, ldwork, v, ldv, &c_b14, c, ldc);
}
/* W := W * V2' */
dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &v[*n - *k], ldv, work, ldwork);
/* C2 := C2 - W */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *m; ++i) {
c[i + (*n - *k + j) * *ldc] -= work[i + j * *ldwork];
}
}
}
}
} else if (lsame_(storev, "R")) {
if (lsame_(direct, "F"))
{
/* Let V = ( V1 V2 ) (V1: first K columns) */
/* where V1 is unit upper triangular. */
if (lsame_(side, "L"))
{
/* Form H * C or H' * C where C = ( C1 ) */
/* ( C2 ) */
/* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
/* W := C1' */
for (j = 0; j < *k; ++j) {
dcopy_(n, &c[j], ldc, &work[j * *ldwork], &c__1);
}
/* W := W * V1' */
dtrmm_("Right", "Upper", "Transpose", "Unit", n, k, &c_b14, v, ldv, work, ldwork);
if (*m > *k)
{
/* W := W + C2'*V2' */
i__1 = *m - *k;
dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &c[*k], ldc, &v[*k * *ldv], ldv, &c_b14, work, ldwork);
}
/* W := W * T' or W * T */
dtrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b14, t, ldt, work, ldwork);
/* C := C - V' * W' */
if (*m > *k) {
/* C2 := C2 - V2' * W' */
i__1 = *m - *k;
dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[*k * *ldv], ldv, work, ldwork, &c_b14, &c[*k], ldc);
}
/* W := W * V1 */
dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14, v, ldv, work, ldwork);
/* C1 := C1 - W' */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *n; ++i) {
c[j + i * *ldc] -= work[i + j * *ldwork];
}
}
} else if (lsame_(side, "R")) {
/* Form C * H or C * H' where C = ( C1 C2 ) */
/* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
/* W := C1 */
for (j = 0; j < *k; ++j) {
dcopy_(m, &c[j * *ldc], &c__1, &work[j * *ldwork], &c__1);
}
/* W := W * V1' */
dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, v, ldv, work, ldwork);
if (*n > *k)
{
/* W := W + C2 * V2' */
i__1 = *n - *k;
dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14,
&c[*k * *ldc], ldc, &v[*k **ldv], ldv, &c_b14, work, ldwork);
}
/* W := W * T or W * T' */
dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, t, ldt, work, ldwork);
/* C := C - W * V */
if (*n > *k)
{
/* C2 := C2 - W * V2 */
i__1 = *n - *k;
dgemm_("No transpose", "No transpose", m, &i__1, k, &c_b25,
work, ldwork, &v[*k * *ldv], ldv, &c_b14, &c[*k * *ldc], ldc);
}
/* W := W * V1 */
dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14, v, ldv, work, ldwork);
/* C1 := C1 - W */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *m; ++i) {
c[i + j * *ldc] -= work[i + j * *ldwork];
}
}
}
} else {
/* Let V = ( V1 V2 ) (V2: last K columns) */
/* where V2 is unit lower triangular. */
if (lsame_(side, "L")) {
/* Form H * C or H' * C where C = ( C1 ) */
/* ( C2 ) */
/* W := C' * V' = (C1'*V1' + C2'*V2') (stored in WORK) */
/* W := C2' */
for (j = 0; j < *k; ++j) {
dcopy_(n, &c[*m - *k + j], ldc, &work[j * *ldwork], &c__1);
}
/* W := W * V2' */
dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &v[(*m - *k) * *ldv], ldv, work, ldwork);
if (*m > *k) {
/* W := W + C1'*V1' */
i__1 = *m - *k;
dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, c, ldc, v, ldv, &c_b14, work, ldwork);
}
/* W := W * T' or W * T */
dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, t, ldt, work, ldwork);
/* C := C - V' * W' */
if (*m > *k) {
/* C1 := C1 - V1' * W' */
i__1 = *m - *k;
dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, v, ldv, work, ldwork, &c_b14, c, ldc);
}
/* W := W * V2 */
dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14, &v[(*m - *k) * *ldv], ldv, work, ldwork);
/* C2 := C2 - W' */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *n; ++i) {
c[*m - *k + j + i * *ldc] -= work[i + j * *ldwork];
}
}
} else if (lsame_(side, "R")) {
/* Form C * H or C * H' where C = ( C1 C2 ) */
/* W := C * V' = (C1*V1' + C2*V2') (stored in WORK) */
/* W := C2 */
for (j = 0; j < *k; ++j) {
dcopy_(m, &c[(*n - *k + j) * *ldc], &c__1, &work[j * *ldwork], &c__1);
}
/* W := W * V2' */
dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &v[(*n - *k) * *ldv], ldv, work, ldwork);
if (*n > *k) {
/* W := W + C1 * V1' */
i__1 = *n - *k;
dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, c, ldc, v, ldv, &c_b14, work, ldwork);
}
/* W := W * T or W * T' */
dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, t, ldt, work, ldwork);
/* C := C - W * V */
if (*n > *k) {
/* C1 := C1 - W * V1 */
i__1 = *n - *k;
dgemm_("No transpose", "No transpose", m, &i__1, k, &c_b25, work, ldwork, v, ldv, &c_b14, c, ldc);
}
/* W := W * V2 */
dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14, &v[(*n - *k) * *ldv], ldv, work, ldwork);
/* C1 := C1 - W */
for (j = 0; j < *k; ++j) {
for (i = 0; i < *m; ++i) {
c[i + (*n - *k + j) * *ldc] -= work[i + j * *ldwork];
}
}
}
}
}
} /* dlarfb_ */
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?