|
|
Ublas : |
From: Gunter Winkler (guwi17_at_[hidden])
Date: 2005-05-24 08:04:26
On Tuesday 24 May 2005 10:47, Andrew Rieck wrote:
> Given a symmetric matrix S and a dense matrix D, what would be the most
> efficient way to construct the matrix M=D^T S D (which will again be
> symmetric) using uBLAS. The former query applies equally to any symmetric
> matrix manipulations which involve non-symmetric matrices and where the
> resulting matrix is symmetric.
>
> Of the three standard techniques that come to mind:
> 1) direct element manipulations, and
This is probably the fastest way, but you have to write the routine (similar
to the ones in operation.hpp)
> 2)
> matrix<double> matrixD(n,n);
> symmetric_matrix<double,lower> matrixS(n,n);
> ...
> symmetric_matrix<double,lower>
> matrixM(prod(trans(matrixD),prod(matrixS,matrixD)));
this prod is not permitted, you have to specify a temporary.
>
> 3) the axpy_prod equivalent of 2) involving a temporary matrix
> matrix<double> matrixTmp(prod(matrixS,matrixD));
> symmetric_matrix<double,lower> matrixM(prod(trans(matrixD),matrixTmp));
AFAIK there is no optimized product to compute this, but you could try this
algorithm:
A = prod(S, D);
R = triangle_adaptor<MATRIX,lower>( prod(trans(D), A )
and R is the lower triangular part of the symmetric result. (Both operations
could use BLAS if available). But I am still not sure if triangle_adaptor
works with any matrix type except dense ...
mfg
Gunter