#include "blaswrap.h" #include "f2c.h" /* Subroutine */ int slarzb_(char *side, char *trans, char *direct, char * storev, integer *m, integer *n, integer *k, integer *l, real *v, integer *ldv, real *t, integer *ldt, real *c__, integer *ldc, real * work, integer *ldwork ) { /* -- LAPACK routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= SLARZB applies a real block reflector H or its transpose H**T to a real distributed M-by-N C from the left or the right. Currently, only STOREV = 'R' and DIRECT = 'B' are supported. Arguments ========= SIDE (input) CHARACTER*1 = 'L': apply H or H' from the Left = 'R': apply H or H' from the Right TRANS (input) CHARACTER*1 = 'N': apply H (No transpose) = 'C': apply H' (Transpose) DIRECT (input) CHARACTER*1 Indicates how H is formed from a product of elementary reflectors = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet) = 'B': H = H(k) . . . H(2) H(1) (Backward) STOREV (input) CHARACTER*1 Indicates how the vectors which define the elementary reflectors are stored: = 'C': Columnwise (not supported yet) = 'R': Rowwise M (input) INTEGER The number of rows of the matrix C. N (input) INTEGER The number of columns of the matrix C. K (input) INTEGER The order of the matrix T (= the number of elementary reflectors whose product defines the block reflector). L (input) INTEGER The number of columns of the matrix V containing the meaningful part of the Householder reflectors. If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. V (input) REAL array, dimension (LDV,NV). If STOREV = 'C', NV = K; if STOREV = 'R', NV = L. LDV (input) INTEGER The leading dimension of the array V. If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K. T (input) REAL array, dimension (LDT,K) The triangular K-by-K matrix T in the representation of the block reflector. LDT (input) INTEGER The leading dimension of the array T. LDT >= K. C (input/output) REAL array, dimension (LDC,N) On entry, the M-by-N matrix C. On exit, C is overwritten by H*C or H'*C or C*H or C*H'. LDC (input) INTEGER The leading dimension of the array C. LDC >= max(1,M). WORK (workspace) REAL array, dimension (LDWORK,K) LDWORK (input) INTEGER The leading dimension of the array WORK. If SIDE = 'L', LDWORK >= max(1,N); if SIDE = 'R', LDWORK >= max(1,M). Further Details =============== Based on contributions by A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA ===================================================================== Quick return if possible Parameter adjustments */ /* Table of constant values */ static integer c__1 = 1; static real c_b13 = 1.f; static real c_b23 = -1.f; /* System generated locals */ integer c_dim1, c_offset, t_dim1, t_offset, v_dim1, v_offset, work_dim1, work_offset, i__1, i__2; /* Local variables */ static integer i__, j, info; extern logical lsame_(char *, char *); extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, integer *, real *, real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *), strmm_(char *, char *, char *, char *, integer *, integer *, real *, real *, integer *, real *, integer *), xerbla_(char *, integer *); static char transt[1]; v_dim1 = *ldv; v_offset = 1 + v_dim1; v -= v_offset; t_dim1 = *ldt; t_offset = 1 + t_dim1; t -= t_offset; c_dim1 = *ldc; c_offset = 1 + c_dim1; c__ -= c_offset; work_dim1 = *ldwork; work_offset = 1 + work_dim1; work -= work_offset; /* Function Body */ if (*m <= 0 || *n <= 0) { return 0; } /* Check for currently supported options */ info = 0; if (! lsame_(direct, "B")) { info = -3; } else if (! lsame_(storev, "R")) { info = -4; } if (info != 0) { i__1 = -info; xerbla_("SLARZB", &i__1); return 0; } if (lsame_(trans, "N")) { *(unsigned char *)transt = 'T'; } else { *(unsigned char *)transt = 'N'; } if (lsame_(side, "L")) { /* Form H * C or H' * C W( 1:n, 1:k ) = C( 1:k, 1:n )' */ i__1 = *k; for (j = 1; j <= i__1; ++j) { scopy_(n, &c__[j + c_dim1], ldc, &work[j * work_dim1 + 1], &c__1); /* L10: */ } /* W( 1:n, 1:k ) = W( 1:n, 1:k ) + ... C( m-l+1:m, 1:n )' * V( 1:k, 1:l )' */ if (*l > 0) { sgemm_("Transpose", "Transpose", n, k, l, &c_b13, &c__[*m - *l + 1 + c_dim1], ldc, &v[v_offset], ldv, &c_b13, &work[ work_offset], ldwork); } /* W( 1:n, 1:k ) = W( 1:n, 1:k ) * T' or W( 1:m, 1:k ) * T */ strmm_("Right", "Lower", transt, "Non-unit", n, k, &c_b13, &t[ t_offset], ldt, &work[work_offset], ldwork); /* C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )' */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *k; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] -= work[j + i__ * work_dim1]; /* L20: */ } /* L30: */ } /* C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ... V( 1:k, 1:l )' * W( 1:n, 1:k )' */ if (*l > 0) { sgemm_("Transpose", "Transpose", l, n, k, &c_b23, &v[v_offset], ldv, &work[work_offset], ldwork, &c_b13, &c__[*m - *l + 1 + c_dim1], ldc); } } else if (lsame_(side, "R")) { /* Form C * H or C * H' W( 1:m, 1:k ) = C( 1:m, 1:k ) */ i__1 = *k; for (j = 1; j <= i__1; ++j) { scopy_(m, &c__[j * c_dim1 + 1], &c__1, &work[j * work_dim1 + 1], & c__1); /* L40: */ } /* W( 1:m, 1:k ) = W( 1:m, 1:k ) + ... C( 1:m, n-l+1:n ) * V( 1:k, 1:l )' */ if (*l > 0) { sgemm_("No transpose", "Transpose", m, k, l, &c_b13, &c__[(*n - * l + 1) * c_dim1 + 1], ldc, &v[v_offset], ldv, &c_b13, & work[work_offset], ldwork); } /* W( 1:m, 1:k ) = W( 1:m, 1:k ) * T or W( 1:m, 1:k ) * T' */ strmm_("Right", "Lower", trans, "Non-unit", m, k, &c_b13, &t[t_offset] , ldt, &work[work_offset], ldwork); /* C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k ) */ i__1 = *k; for (j = 1; j <= i__1; ++j) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { c__[i__ + j * c_dim1] -= work[i__ + j * work_dim1]; /* L50: */ } /* L60: */ } /* C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ... W( 1:m, 1:k ) * V( 1:k, 1:l ) */ if (*l > 0) { sgemm_("No transpose", "No transpose", m, l, k, &c_b23, &work[ work_offset], ldwork, &v[v_offset], ldv, &c_b13, &c__[(*n - *l + 1) * c_dim1 + 1], ldc); } } return 0; /* End of SLARZB */ } /* slarzb_ */