[eigen] Re: efficient trace of product of symmetric matrices

[Patrick <pat05.mlst@xxxxxxxxx>]

Hi everyone,

I'm wondering whether it is currently possible to compute

tr(A*B),   where   A \in Sym(n,n), B \in Sym(n,n)

efficiently using the development branch of Eigen3. 

Since A and B are symmetric and 

tr(A*B) = sum_ij a_ij * b_ij,

the best way I can think of doing this is 

double tr=0;

for(size_t i=0; i<A.rows(); ++i)

{

    tr += A(i,i) * B(i,i);

    for(size_t j=0; j<i; ++j)

        tr += 2 * A(i,j)*B(i,j);

}

I noticed that the triangularView methods do not expose the array 

interface, so am I correct in assuming I cannot vectorise this using 

triangular views?

Furthermore, when trying to make a matrix symmetric by calling 

A.triangularView<StrictlyUpper>() = A.triangularView<StrictlyLower>(),

I noticed that no assignment was taking place, i.e. the matrix A remained as is.

Is this the expected behaviour?

Thanks.

Kind regards,

Patrick