Re: [eigen] A not so simple product |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] A not so simple product*From*: "Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>*Date*: Tue, 9 Dec 2008 15:32:39 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:received:message-id:date:from:to :subject:in-reply-to:mime-version:content-type :content-transfer-encoding:content-disposition:references; bh=hrC5wwitZFfVKBFOBUMOgZVmxf3/AV6dF63A20pD3NQ=; b=cYaqzmjMfDyvsZgDEldY/qClYX9z2qX0CfCqzhelJrx0fJtKfZeoKtz5vjdsGRGwwu 6VcO3/SvZI7hxpSZwkMouZkV4KYuhZxM3h4Zdw0nldWK8Wpl5/KaJGsb2/ww19Qb5fkB J3kfJGLoEHbetA40odP5eij4OOtEd6ACp9cUE=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=message-id:date:from:to:subject:in-reply-to:mime-version :content-type:content-transfer-encoding:content-disposition :references; b=DNQYICeAcNkfxvg5ZjuIkbIwZdX613IFViEPN5tMKpoESnVeoYguEImbCurrNNetdk kgg+bAA0xmR9kX2Tpju5Lk9UwyEIuglfohMg0BvRpheAwcb9LZaXyxFWjrMyA/jBsOdi 6671aIRJXLayG86PKnWMRRV5D2+84nJaSebdQ=

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 >>> >>> --- >>> >> >> > ---

**Follow-Ups**:**Re: [eigen] A not so simple product***From:*Benjamin Schindler

**References**:**[eigen] A not so simple product***From:*Benjamin Schindler

**Re: [eigen] A not so simple product***From:*Gael Guennebaud

**Re: [eigen] A not so simple product***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] A not so simple product** - Next by Date:
**Re: [eigen] A not so simple product** - Previous by thread:
**Re: [eigen] A not so simple product** - Next by thread:
**Re: [eigen] A not so simple product**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |