zlarfb.c

来自「InsightToolkit-1.4.0(有大量的优化算法程序)」· C语言 代码 · 共 587 行 · 第 1/2 页

                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i += work[i__2].i;
                    }
                }

            } 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) {
                    zcopy_(m, &c[(*n - *k + j) * *ldc], &c__1, &work[j* *ldwork], &c__1);
                }

/*              W := W * V2 */

                ztrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b15, &v[*n - *k], ldv, work, ldwork);
                if (*n > *k) {

/*                 W := W + C1 * V1 */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "No transpose", m, k, &i__1, &c_b15,
                           c, ldc, v, ldv, &c_b15, work, ldwork);
                }

/*              W := W * T  or  W * T' */

                ztrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b15, t, ldt, work, ldwork);

/*              C := C - W * V' */

                if (*n > *k) {

/*                 C1 := C1 - W * V1' */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "Conjugate transpose", m, &i__1, k,
                           &c_b26, work, ldwork, v, ldv, &c_b15, c, ldc);
                }

/*              W := W * V2' */

                ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", m, k, &c_b15, &v[*n - *k], ldv, work, ldwork);

/*              C2 := C2 - W */

                for (j = 0; j < *k; ++j) {
                    for (i = 0; i < *m; ++i) {
                        i__1 = i + (*n - *k + j) * *ldc;
                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i -= work[i__2].i;
                    }
                }
            }
        }

    } 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) {
                    zcopy_(n, &c[j], ldc, &work[j* *ldwork], &c__1);
                    zlacgv_(n, &work[j* *ldwork], &c__1);
                }

/*              W := W * V1' */

                ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", n, k, &c_b15, v, ldv, work, ldwork);
                if (*m > *k) {

/*                 W := W + C2'*V2' */

                    i__1 = *m - *k;
                    zgemm_("Conjugate transpose", "Conjugate transpose", n, k,
                           &i__1, &c_b15, &c[*k], ldc, &v[*k * *ldv], ldv, &c_b15, work, ldwork);
                }

/*              W := W * T'  or  W * T */

                ztrmm_("Right", "Upper", transt, "Non-unit", n, k, &c_b15, t, ldt, work, ldwork);

/*              C := C - V' * W' */

                if (*m > *k) {

/*                 C2 := C2 - V2' * W' */

                    i__1 = *m - *k;
                    zgemm_("Conjugate transpose", "Conjugate transpose", &i__1,
                           n, k, &c_b26, &v[*k * *ldv], ldv, work, ldwork, &c_b15, &c[*k], ldc);
                }

/*              W := W * V1 */

                ztrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b15, v, ldv, work, ldwork);

/*              C1 := C1 - W' */

                for (j = 0; j < *k; ++j) {
                    for (i = 0; i < *n; ++i) {
                        i__1 = j + i * *ldc;
                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i += work[i__2].i;
                    }
                }

            } 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) {
                    zcopy_(m, &c[j * *ldc], &c__1, &work[j* *ldwork], &c__1);
                }

/*              W := W * V1' */

                ztrmm_("Right", "Upper", "Conjugate transpose", "Unit", m, k, &c_b15, v, ldv, work, ldwork);
                if (*n > *k) {

/*                 W := W + C2 * V2' */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "Conjugate transpose", m, k, &i__1,
                           &c_b15, &c[*k * *ldc], ldc, &v[*k * *ldv], ldv, &c_b15, work, ldwork);
                }

/*              W := W * T  or  W * T' */

                ztrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b15, t, ldt, work, ldwork);

/*              C := C - W * V */

                if (*n > *k) {

/*                 C2 := C2 - W * V2 */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "No transpose", m, &i__1, k, &c_b26,
                           work, ldwork, &v[*k * *ldv], ldv, &c_b15, &c[*k * *ldc], ldc);
                }

/*              W := W * V1 */

                ztrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b15, v, ldv, work, ldwork);

/*              C1 := C1 - W */

                for (j = 0; j < *k; ++j) {
                    for (i = 0; i < *m; ++i) {
                        i__1 = i + j * *ldc;
                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i -= work[i__2].i;
                    }
                }
            }

        } 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) {
                    zcopy_(n, &c[*m - *k + j], ldc, &work[j* *ldwork], &c__1);
                    zlacgv_(n, &work[j* *ldwork], &c__1);
                }

/*              W := W * V2' */

                ztrmm_("Right", "Lower", "Conjugate transpose", "Unit", n, k, &c_b15, &v[(*m - *k) * *ldv], ldv, work, ldwork);
                if (*m > *k) {

/*                 W := W + C1'*V1' */

                    i__1 = *m - *k;
                    zgemm_("Conjugate transpose", "Conjugate transpose", n, k,
                           &i__1, &c_b15, c, ldc, v, ldv, &c_b15, work, ldwork);
                }

/*              W := W * T'  or  W * T */

                ztrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b15, t, ldt, work, ldwork);

/*              C := C - V' * W' */

                if (*m > *k) {

/*                 C1 := C1 - V1' * W' */

                    i__1 = *m - *k;
                    zgemm_("Conjugate transpose", "Conjugate transpose",
                           &i__1, n, k, &c_b26, v, ldv, work, ldwork, &c_b15, c, ldc);
                }

/*              W := W * V2 */

                ztrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b15, &v[(*m - *k) * *ldv], ldv, work, ldwork);

/*              C2 := C2 - W' */

                for (j = 0; j < *k; ++j) {
                    for (i = 0; i < *n; ++i) {
                        i__1 = *m - *k + j + i * *ldc;
                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i += work[i__2].i;
                    }
                }

            } 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) {
                    zcopy_(m, &c[(*n - *k + j) * *ldc], &c__1, &work[j* *ldwork], &c__1);
                }

/*              W := W * V2' */

                ztrmm_("Right", "Lower", "Conjugate transpose", "Unit", m, k, &c_b15, &v[(*n - *k) * *ldv], ldv, work, ldwork);
                if (*n > *k) {

/*                 W := W + C1 * V1' */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "Conjugate transpose", m, k, &i__1,
                           &c_b15, c, ldc, v, ldv, &c_b15, work, ldwork);
                }

/*              W := W * T  or  W * T' */

                ztrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b15, t, ldt, work, ldwork);

/*              C := C - W * V */

                if (*n > *k) {

/*                 C1 := C1 - W * V1 */

                    i__1 = *n - *k;
                    zgemm_("No transpose", "No transpose", m, &i__1, k,
                           &c_b26, work, ldwork, v, ldv, &c_b15, c, ldc);
                }

/*              W := W * V2 */

                ztrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b15, &v[(*n - *k) * *ldv], ldv, work, ldwork);

/*              C1 := C1 - W */

                for (j = 0; j < *k; ++j) {
                    for (i = 0; i < *m; ++i) {
                        i__1 = i + (*n - *k + j) * *ldc;
                        i__2 = i + j * *ldwork;
                        c[i__1].r -= work[i__2].r, c[i__1].i -= work[i__2].i;
                    }
                }
            }
        }
    }
} /* zlarfb_ */