BLAS operation 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.

• Is this part of an iterative process or done only once? If part of lop, does U or B or both change in each iteration? Jun 16 '19 at 1:26
• I added more detail in the post on what I'm trying to do
– vibe
Jun 16 '19 at 2:17

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.

• thanks, I've considered this but I'd prefer to avoid allocating more workspace if possible
– vibe
Jun 16 '19 at 2:20