dlarfb.c

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

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

/*              W := W * V1' */

/*<    >*/
                dtrmm_("Right", "Lower", "Transpose", "Unit", m, k, &c_b14, &
                        v[v_offset], ldv, &work[work_offset], ldwork, (ftnlen)
                        5, (ftnlen)5, (ftnlen)9, (ftnlen)4);

/*              C1 := C1 - W */

/*<                DO 60 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 50 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];
/*<    50             CONTINUE >*/
/* L50: */
                    }
/*<    60          CONTINUE >*/
/* L60: */
                }
/*<             END IF >*/
            }

/*<          ELSE >*/
        } else {

/*           Let  V =  ( V1 ) */
/*                     ( V2 )    (last K rows) */
/*           where  V2  is unit upper 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 70 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);
/*<    70          CONTINUE >*/
/* L70: */
                }

/*              W := W * V2 */

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

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

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

/*              W := W * V2' */

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

/*              C2 := C2 - W' */

/*<                DO 90 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 80 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];
/*<    80             CONTINUE >*/
/* L80: */
                    }
/*<    90          CONTINUE >*/
/* L90: */
                }

/*<             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 100 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);
/*<   100          CONTINUE >*/
/* L100: */
                }

/*              W := W * V2 */

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

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

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

/*              W := W * V2' */

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

/*              C2 := C2 - W */

/*<                DO 120 J = 1, K >*/
                i__1 = *k;
                for (j = 1; j <= i__1; ++j) {
/*<                   DO 110 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];
/*<   110             CONTINUE >*/
/* L110: */
                    }
/*<   120          CONTINUE >*/
/* L120: */
                }
/*<             END IF >*/
            }
/*<          END IF >*/
        }

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

/*<          IF( LSAME( DIRECT, 'F' ) ) THEN >*/
        if (lsame_(direct, "F", (ftnlen)1, (ftnlen)1)) {

/*           Let  V =  ( V1  V2 )    (V1: first K columns) */
/*           where  V1  is unit upper 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 := C1' */

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

/*              W := W * V1' */

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

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

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

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

/*<    >*/
                dtrmm_("Right", "Upper", 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) {

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

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