1

Given two matrices, A and B, where B is symetric (and positive semi-definite), What is the best (fastest) way to calculate A`*B*A?

Currently, using BLAS, I first compute C=B*A using dsymm (introducing a temporary matrix C) and then A`*C using dgemm.

Is there a better (faster, no temporaries) way to do this using BLAS and mkl?

Thanks.

nir
  • 96
  • 3
  • 1
    See also http://stackoverflow.com/questions/11139933/efficient-way-of-computing-matrix-product-axa and http://stackoverflow.com/questions/13292325/m-x-s-x-tm-multiplication-with-blas – finnw Feb 26 '13 at 18:50
  • Thanks @finnw. However, there is no answer in the posts you mention. – nir Feb 27 '13 at 05:07
  • 1
    This is the third time I see this question with slight modifications asked in SO (first one was my own question), and the lack of answers have led me to believe that your way own suggestion is fastest possible. – Jouni Helske Feb 28 '13 at 17:24

1 Answers1

1

I'll offer somekind of answer: Compared to the general case A*B*C you know that the end result is symmetric matrix. After computing C=B*A with BLAS subroutine dsymm, you want to compute A'C, but you only need to compute the upper diagonal part of the matrix and the copy the strictly upper diagonal part to the lower diagonal part.

Unfortunately there doesn't seem to be a BLAS routine where you can claim beforehand that given two general matrices, the output matrix will be symmetric. I'm not sure if it would be beneficial to write you own function for this. This probably depends on the size of your matrices and the implementation.

EDIT: This idea seems to be addressed recently here: A Matrix Multiplication Routine that Updates Only the Upper or Lower Triangular Part of the Result Matrix

Jouni Helske
  • 6,427
  • 29
  • 52
  • This link indeed seems interesting. Alas, it only saves about 20% of the runtime, and in my specific application I still need to copy the results to the lower triangle - which leaves me with almost no gain... Thanks. – nir Mar 03 '13 at 09:13