2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Brian Borchers Jun 16 at 1:26
  • $\begingroup$ I added more detail in the post on what I'm trying to do $\endgroup$ – vibe Jun 16 at 2:17
2
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.