[eigen] SelfAdjointView<> as LHS expression for optimizations?

[Jose Luis Blanco <joseluisblancoc@xxxxxxxxx>]

Hi all,

The operation I have in mind is:
  RESULT = H * C * H.transpose()
with C symmetric and squared.

The point is that two "problems" (well... I have to admit, small
ones!) occur with a naive Eigen expresion "H * C * H.transpose()":

1) Many multiplications could be saved by realizing that the result is
exactly symmetric. I know the implications may be far reaching, but,
could SelfAdjointView<> be extended so it can be used to store
results? I mean, something like:

Eigen::MatrixXd RESULT;
RESULT.selfadjointView<Eigen::Lower>() =
H*C.selfadjointView<Eigen::Lower>()*H.transpose();

or an alternative that wouldn't require LHS  = operator for SelfAdjointView<>:

Eigen::MatrixXd RESULT =
(H*C.selfadjointView<Eigen::Lower>()*H.transpose()).selfAdjointView<Eigen::Lower>();

Any ideas about the difficulty of implementing any of these two ideas?
Or... is there anything similar already done?

2) Apart from efficiency, there is this numerical accuracy problem of
the resulting matrix not being *exactly* symmetric... a*b*c != c*b*a
for most architectures. Ok, the errors are tiny, but it may lead to
inconsistencies that some time ago affected the Cholesky
implementation I used to use... that's why I worry about all this in
the first place!


Best,
JL