Re: [boost] [gsoc18][ublas] Proposal to add advanced matrix operations

Subject: Re: [boost] [gsoc18][ublas] Proposal to add advanced matrix operations
From: Hans Dembinski (hans.dembinski_at_[hidden])
Date: 2018-01-23 13:26:24

> On 22. Jan 2018, at 23:07, Marcin Copik <mcopik_at_[hidden]> wrote:
> 
> 
> 
> pon., 22 sty 2018 o 11:16 użytkownik Hans Dembinski via Boost <boost_at_[hidden] <mailto:boost_at_[hidden]>> napisał:
> Hi David, Rajaditya, Artyom, Shikhar,
> 
> > On 21. Jan 2018, at 13:00, David Bellot via Boost <boost_at_[hidden] <mailto:boost_at_[hidden]>> wrote:
> >
> 
> I used the library in projects and it is very convenient. You just write your expressions as you would on paper and let Eigen generate the fastest possible code out of that. I only dislike that Eigen does not play well with "auto" since C++11, because they use the assignment to a matrix or vector to trigger the evaluation of the expression template.
> 
> As far as I know, direct usage of auto is not possible in any Expression Templates library since the assigment triggers an evaluation, as you already noticed. It's a minor inconviencen since it requires from the you user to know if final type of an expression should be a matrix or a vector. It might be a larger problem in a library which provides specialized types for matrices of different shapes, such as: symmetric, upper/lower triangular, diagonal. Such types are benefitial since it helps the library to find a more optimized kernel, e.g. use trmm instead of gemm for a matrix-matrix product if one of operands is lower/upper triangular.
> Then, the return type of an expression depends on types of input matrices and the actual computation and one should not expect from the user to be able to guess it correctly. As far as I know, only Blaze implements specialized types called adaptors (Eigen has views but their applicability is rather limited).

That is a good point. In Eigen, you can use .eval(), see https://eigen.tuxfamily.org/dox/TopicPitfalls.html

Example:
MatrixXd A, B;
auto C_bad = A * B; // this is an expression template, not to be used
auto C_good = (A * B).eval(); // this is a matrix type, eval() can chose the best type to represent the result

> Fortunately, there should be an easy solution available in each library. Each ET library should have a return type implemented in each expression node. Then, it is possible to implement a simple evaluate function:
> template<typename Expr>
> typename Expr::expression_type evaluate(Expr && expr)
> {
>     return typename Expr::expression_type(expr);
> }
> It should always enforce the evaluation of an expression and one can safely return an expression from a function without knowing the actual type.

Agreed. I am not sure whether a free function or a method call is better. Apparently, the Eigen guys went for the latter.

Best regards,
Hans