dlarfb.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*              W := W * V1 */

/*<    >*/
                dtrmm_("Right", "Upper", "No transpose", "Unit", n, k, &c_b14,
                         &v[v_offset], ldv, &work[work_offset], ldwork, (
                        ftnlen)5, (ftnlen)5, (ftnlen)12, (ftnlen)4);

/*              C1 := C1 - W' */

/*<                DO 150 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 140 I = 1, N >*/
                    i__2 = *n;
                    for (i__ = 1; i__ <= i__2; ++i__) {
/*<                      C( J, I ) = C( J, I ) - WORK( I, J ) >*/
                        c__[j + i__ * c_dim1] -= work[i__ + j * work_dim1];
/*<   140             CONTINUE >*/
/* L140: */
                    }
/*<   150          CONTINUE >*/
/* L150: */
                }

/*<             ELSE IF( LSAME( SIDE, 'R' ) ) THEN >*/
            } else if (lsame_(side, "R", (ftnlen)1, (ftnlen)1)) {

/*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */

/*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */

/*              W := C1 */

/*<                DO 160 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 ) >*/
                    dcopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j * 
                            work_dim1 + 1], &c__1);
/*<   160          CONTINUE >*/
/* L160: */
                }

/*              W := W * V1' */

/*<    >*/
                dtrmm_("Right", "Upper", "Transpose", "Unit", m, k, &c_b14, &
                        v[v_offset], ldv, &work[work_offset], ldwork, (ftnlen)
                        5, (ftnlen)5, (ftnlen)9, (ftnlen)4);
/*<                IF( N.GT.K ) THEN >*/
                if (*n > *k) {

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

/*<    >*/
                    i__1 = *n - *k;
                    dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
                            c__[(*k + 1) * c_dim1 + 1], ldc, &v[(*k + 1) * 
                            v_dim1 + 1], ldv, &c_b14, &work[work_offset], 
                            ldwork, (ftnlen)12, (ftnlen)9);
/*<                END IF >*/
                }

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

/*<    >*/
                dtrmm_("Right", "Upper", trans, "Non-unit", m, k, &c_b14, &t[
                        t_offset], ldt, &work[work_offset], ldwork, (ftnlen)5,
                         (ftnlen)5, (ftnlen)1, (ftnlen)8);

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

/*<                IF( N.GT.K ) THEN >*/
                if (*n > *k) {

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

/*<    >*/
                    i__1 = *n - *k;
                    dgemm_("No transpose", "No transpose", m, &i__1, k, &
                            c_b25, &work[work_offset], ldwork, &v[(*k + 1) * 
                            v_dim1 + 1], ldv, &c_b14, &c__[(*k + 1) * c_dim1 
                            + 1], ldc, (ftnlen)12, (ftnlen)12);
/*<                END IF >*/
                }

/*              W := W * V1 */

/*<    >*/
                dtrmm_("Right", "Upper", "No transpose", "Unit", m, k, &c_b14,
                         &v[v_offset], ldv, &work[work_offset], ldwork, (
                        ftnlen)5, (ftnlen)5, (ftnlen)12, (ftnlen)4);

/*              C1 := C1 - W */

/*<                DO 180 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 170 I = 1, M >*/
                    i__2 = *m;
                    for (i__ = 1; i__ <= i__2; ++i__) {
/*<                      C( I, J ) = C( I, J ) - WORK( I, J ) >*/
                        c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1];
/*<   170             CONTINUE >*/
/* L170: */
                    }
/*<   180          CONTINUE >*/
/* L180: */
                }

/*<             END IF >*/
            }

/*<          ELSE >*/
        } else {

/*           Let  V =  ( V1  V2 )    (V2: last K columns) */
/*           where  V2  is unit lower triangular. */

/*<             IF( LSAME( SIDE, 'L' ) ) THEN >*/
            if (lsame_(side, "L", (ftnlen)1, (ftnlen)1)) {

/*              Form  H * C  or  H' * C  where  C = ( C1 ) */
/*                                                  ( C2 ) */

/*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK) */

/*              W := C2' */

/*<                DO 190 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   CALL DCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 ) >*/
                    dcopy_(n, &c__[*m - *k + j + c_dim1], ldc, &work[j * 
                            work_dim1 + 1], &c__1);
/*<   190          CONTINUE >*/
/* L190: */
                }

/*              W := W * V2' */

/*<    >*/
                dtrmm_("Right", "Lower", "Transpose", "Unit", n, k, &c_b14, &
                        v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
                        , ldwork, (ftnlen)5, (ftnlen)5, (ftnlen)9, (ftnlen)4);
/*<                IF( M.GT.K ) THEN >*/
                if (*m > *k) {

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

/*<    >*/
                    i__1 = *m - *k;
                    dgemm_("Transpose", "Transpose", n, k, &i__1, &c_b14, &
                            c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
                            work[work_offset], ldwork, (ftnlen)9, (ftnlen)9);
/*<                END IF >*/
                }

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

/*<    >*/
                dtrmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b14, &t[
                        t_offset], ldt, &work[work_offset], ldwork, (ftnlen)5,
                         (ftnlen)5, (ftnlen)1, (ftnlen)8);

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

/*<                IF( M.GT.K ) THEN >*/
                if (*m > *k) {

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

/*<    >*/
                    i__1 = *m - *k;
                    dgemm_("Transpose", "Transpose", &i__1, n, k, &c_b25, &v[
                            v_offset], ldv, &work[work_offset], ldwork, &
                            c_b14, &c__[c_offset], ldc, (ftnlen)9, (ftnlen)9);
/*<                END IF >*/
                }

/*              W := W * V2 */

/*<    >*/
                dtrmm_("Right", "Lower", "No transpose", "Unit", n, k, &c_b14,
                         &v[(*m - *k + 1) * v_dim1 + 1], ldv, &work[
                        work_offset], ldwork, (ftnlen)5, (ftnlen)5, (ftnlen)
                        12, (ftnlen)4);

/*              C2 := C2 - W' */

/*<                DO 210 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 200 I = 1, N >*/
                    i__2 = *n;
                    for (i__ = 1; i__ <= i__2; ++i__) {
/*<                      C( M-K+J, I ) = C( M-K+J, I ) - WORK( I, J ) >*/
                        c__[*m - *k + j + i__ * c_dim1] -= work[i__ + j * 
                                work_dim1];
/*<   200             CONTINUE >*/
/* L200: */
                    }
/*<   210          CONTINUE >*/
/* L210: */
                }

/*<             ELSE IF( LSAME( SIDE, 'R' ) ) THEN >*/
            } else if (lsame_(side, "R", (ftnlen)1, (ftnlen)1)) {

/*              Form  C * H  or  C * H'  where  C = ( C1  C2 ) */

/*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK) */

/*              W := C2 */

/*<                DO 220 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   CALL DCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) >*/
                    dcopy_(m, &c__[(*n - *k + j) * c_dim1 + 1], &c__1, &work[
                            j * work_dim1 + 1], &c__1);
/*<   220          CONTINUE >*/
/* L220: */
                }

/*              W := W * V2' */

/*<    >*/
                dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
                        v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[work_offset]
                        , ldwork, (ftnlen)5, (ftnlen)5, (ftnlen)9, (ftnlen)4);
/*<                IF( N.GT.K ) THEN >*/
                if (*n > *k) {

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

/*<    >*/
                    i__1 = *n - *k;
                    dgemm_("No transpose", "Transpose", m, k, &i__1, &c_b14, &
                            c__[c_offset], ldc, &v[v_offset], ldv, &c_b14, &
                            work[work_offset], ldwork, (ftnlen)12, (ftnlen)9);
/*<                END IF >*/
                }

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

/*<    >*/
                dtrmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b14, &t[
                        t_offset], ldt, &work[work_offset], ldwork, (ftnlen)5,
                         (ftnlen)5, (ftnlen)1, (ftnlen)8);

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

/*<                IF( N.GT.K ) THEN >*/
                if (*n > *k) {

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

/*<    >*/
                    i__1 = *n - *k;
                    dgemm_("No transpose", "No transpose", m, &i__1, k, &
                            c_b25, &work[work_offset], ldwork, &v[v_offset], 
                            ldv, &c_b14, &c__[c_offset], ldc, (ftnlen)12, (
                            ftnlen)12);
/*<                END IF >*/
                }

/*              W := W * V2 */

/*<    >*/
                dtrmm_("Right", "Lower", "No transpose", "Unit", m, k, &c_b14,
                         &v[(*n - *k + 1) * v_dim1 + 1], ldv, &work[
                        work_offset], ldwork, (ftnlen)5, (ftnlen)5, (ftnlen)
                        12, (ftnlen)4);

/*              C1 := C1 - W */

/*<                DO 240 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 230 I = 1, M >*/
                    i__2 = *m;
                    for (i__ = 1; i__ <= i__2; ++i__) {
/*<                      C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J ) >*/
                        c__[i__ + (*n - *k + j) * c_dim1] -= work[i__ + j * 
                                work_dim1];
/*<   230             CONTINUE >*/
/* L230: */
                    }
/*<   240          CONTINUE >*/
/* L240: */
                }

/*<             END IF >*/
            }

/*<          END IF >*/
        }
/*<       END IF >*/
    }

/*<       RETURN >*/
    return 0;

/*     End of DLARFB */

/*<       END >*/
} /* dlarfb_ */

#ifdef __cplusplus
        }
#endif