Re: [eigen] A not so simple product

["Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>]

OK, just had a chat with Gael over IRC.

You should use the diagonal matrix product approach. It will be easy
to add vectorization support for that, I add that to the TODO.

Cheers,
Benoit

2008/12/9 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> Indeed, the "obvious" way
>
> for(i = 0; i < N; i++) m.row(i) /= rho(i)*w(i);
>
> doesn't vectorize because rows have size 3... so if you want
> vectorization you need a complicated trick as Gael shows.
>
> Gael: in order to achieve vectorization of diagonal products, wouldn't
> it be simplest to just change DiagonalProduct.h so it takes in account
> the fact that coefficients come from a vector?
>
> Something like replacing m_lhs.coeff(row, row) by
> m_lhs.m_coeffs.coeff(row), since we know that m_lhs is a
> DiagonalMatrix?
>
> Cheers,
> Benoit
>
> 2008/12/9 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>
>> Hi,
>>
>> I guess there is no M and you wanted to write "VectorXf w(N);" instead of
>> "VectorXf w(M);" ?
>>
>> Anyway, you can write it as a single expression like that:
>>
>> m = ( (rho.cwise()*w).cwise().inverse().asDiagonal() ) * m;
>>
>> the matrix-matrix product will be optimized since Eigen knows that the
>> left-hand side (lhs) matrix is a diagonal matrix. Unfortunately, Eigen does
>> not take into account yet the fact that the diagonal of lhs is actually
>> stored as a vector, and therefore this expression won't be vectorized.
>>
>> The good think is that this is an optimization we can quite easily add. So I
>> would recommand to do so and wait a bit that such expressions get vectorized
>> by Eigen.
>>
>> Note that in the meantime, if performance matters you can also do:
>>
>> rho (or w, or a temporary) = rho.cwise()*w).cwise().inverse();
>> for (int k=0: k<3; ++k)
>>   m.col(i).cwise() *= rho (or w, or a temporary);
>>
>> in which case all computations will be vectorized at the cost of additionnal
>> read/writes operations.
>>
>>
>>
>> Below are some proposals related to that problem:
>>
>>
>> Another option would be to add a mechanism to extend a vector to a matrix
>> with repetitions:
>>
>> m = m .cwise()* ( (rho.cwise()*w).cwise().inverse().duplicate(3) );
>>
>> I'm not sure about "duplicate" for the name but I know matlab offers such a
>> feature. For the vectorization this is a bit tricky because such a
>> "Duplicate" expression can be vectorized in both direction: you can either
>> return column packets of row packets (by repeating a scalar in each
>> component of the packet) to accomodate the storage of m. However, this fine
>> tunning of the packet generation is something we cannot support without big
>> changes in Eigen (or ugly workarounds).
>>
>>
>> Finally, a third solution to implement this kind of thing could be to extend
>> the rowwise/colwise API:
>>
>> m = m.colwise() * (rho.cwise()*w).cwise().inverse();
>>
>> perhaps this last proposal is a bit confusing.
>>
>>
>> In any case the choice for the traversal of m is a bit tricky. Assuming m is
>> column major, we can process it column per column that is good wrt to the
>> cache of m but which requires a temporary to store
>> "rho.cwise()*w).cwise().inverse()", or we process m per block of packet-size
>> rows...
>>
>>
>> cheers,
>> Gael.
>>
>>
>> On Tue, Dec 9, 2008 at 2:10 PM, Benjamin Schindler
>> <bschindler@xxxxxxxxxxxxxxx> wrote:
>>>
>>> Hi
>>>
>>> Okay, I'm trying to figure the most efficient way to vectorize this
>>> product:
>>>
>>> MatrixXf m(N,3);
>>> VectorXf rho(N);
>>> VectorXf w(M);
>>>
>>> What I would like to do: The i'th row of m is suppposed to get multiplied
>>> by 1 / (rho(i)*w(i))
>>>
>>> I've not yet found a way which actually works as a short
>>> eigen-expression... is there any or do I have to go back to the old rusty
>>> for loop? :)
>>>
>>> Thanks
>>> Benjamin
>>>
>>> ---
>>>
>>
>>
>

---