mkl/blas routine for C=AA'B

[may_ka]

Hi there,

is there any mkl/blas function which performs the operation C=AA'B in on go. Currently I use an intermediate array T and dgemm: T=A'B;C=AT'. I am wondering whether there is a more efficient way since A is always the same matrix.

Thanks.

[mecej4]

Did you mean to write C = AT rather than C = AT' ?

If A is "always the same", you could form T = AA' just once, using ?GEMM. Then, do the ?GEMM operation C = TB for each B as you go.

A more useful answer could be given if you describe how C is going to be used. Are the matrices square?

[may_ka]

Hi,

Thanks for the response.

I may have mixed the dimensions/transpositions .............. but in essence T is such that it can serve as an intermediate matrix.

A is always the same. Generating AA' would simply be a rank update (?syrk) yielding a symmetric matrix. However, AA' is not feasible because of its dimensions. A can be anything but can have row dimension of several hundreds of thousands. In contrast its column dimension may be only several 10th of thousand. The same holds for B.

The whole system is accutally a nice example for a situation where T=A'B and subsequently C=AT' is much faster (even when done repeatedly) than forming T=AA' once (if feasible) and doing C=TB.

C must have the same dimension as matrix B. It is neither squared nor symmetric.

My thoughts were since A is the same in both operations and C must be the same dimension as B, there might be a more efficient way than putting A through the CPU twice.

Cheers