BLAS operation question

Question

I want to perform the following operation: $$ A = A + U B^T $$ where $A$ is $m \times n$ dense, $U$ is $m \times m$ upper triangular, and $B$ is $n \times m$ dense. The BLAS function dtrmm comes close but cannot handle $B^T$ and also does not compute the sum part. Does anyone see a clever way to compute this with Level 3 BLAS calls?

EDIT: some more details on what I'm doing.

What I am doing is computing the quantity $U L^T$ in-place, where $U$ is upper triangular, and $L$ is unit lower triangular. I have a $N \times N$ matrix $A$ where $U$ is stored in the upper triangle (plus diagonal), and $L$ is stored in the lower triangle. I want to store the product $U L^T$ in-place in the upper triangle of the matrix $A$. I am using a recursive algorithm which works as follows: $$ U L^T = \pmatrix{U_{11} & U_{12} \\ 0 & U_{22}} \pmatrix{L_{11}^T & L_{21}^T \\ 0 & L_{22}^T} = \pmatrix{U_{11} L_{11}^T & U_{11} L_{21}^T + U_{12} L_{22}^T \\ 0 & U_{22} L_{22}^T} $$ So my algorithm proceeds as follows:

1. TRMM: $U_{12} := U_{12} L_{22}^T$
2. ????: $U_{12} := U_{12} + U_{11} L_{21}^T$ (this is my question above)
3. recursion: call routine recursively to compute $U_{11} L_{11}^T$ in-place
4. recursion: call routine recursively to compute $U_{22} L_{22}^T$ in-place

So, step 2 is where I can't find an efficient method. Allocating a temporary workspace matrix is possible, as noted in Anton's answer, but its not ideal since many allocations would be required due to the recursion.

EDIT2:

One possibility I just realized is:

2(a): TRMM: $L_{21} := L_{21} U_{11}^T$

2(b): ADD: $U_{12} := U_{12} + L_{21}^T$

This would destroy the matrix $L$ but I think it would work, and not require additional workspace to be allocated.

Answer 1

Assuming we are looking at standard BLAS interface, I think you cannot really take advantage of the upper-triangular structure of matrix $U$, so you are back to a general matrix-matrix product via ?gemm.

You can try sticking to the triangular matrix-matrix product, but I doubt it is going to be faster.

1. Create a temporary matrix $C=B^T$, using, say, a for loop and a ?copy functionality with strides (for transposition).
2. Use ?trmm to perform $C=UC$
3. Update $A$ using ?axpy to perform $A=A+C$
4. Destroy a temporary matrix $C$

This approach requires to allocate additional memory for matrix $C\in\mathbb F^{m\times n}$, and perform $(n+1)$ calls to level 1 BLAS subroutines ($n$ ?copys and 1 ?axpy) in addition to a level 3 BLAS triangular matrix-matrix product.

Note, some BLAS/LAPACK implementations have their extensions (or even variations on the standard function). For example, see cublas?trmm which comes much closer to what you want.