Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower> |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>*From*: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>*Date*: Tue, 15 Jun 2010 16:29:12 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type :content-transfer-encoding; bh=R4rW5awTcCqLAnDPpBOuCb4AWKbJ0L/pDkZx8q8lxRQ=; b=xcvSXEv3ZxcFrmmjVmqudZsv1Qy4XF0oL0PcdsL1tukcvB/aHfM3ldcFsbUwfYxWvn FjDxbyquZ7h/rTcVboHyGuIBrg924PTyAil3/jGEtfc05BWtcvSxEBiEXT75oeVsXm4t SC3pxMMD4/8u5/M4vVRi8eBxk0NYCGZieuE94=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; b=c42x/xogXI3V4hAKGuM4zFP+zw/ebCsF5N54ijE/UL+c386sx1T1xbVwggTiTtVOJV ciy0YdsIcuMYvsIA4/Pf2wnn9oJD9ZXrIBR48cUQOnWcq/+9MetRl2KP3MScVBjbPD3G YzUBKniJlbkqywYJJwWXUA+BIx7zaO6P+ovS4=

Actually, the pb is line 50 where you store a reference to the temporary matrix returned by BlockHouseholderTriangularFactor. I would rather write it like this: const TriangularView<Derived,UnitLower>& Vtri(V); TriangularFactorType T = BlockHouseholderTriangularFactor(Vtmp, coeffs,V); A -= Vtri * T.transpose() * Vtri.transpose() * A; But here there are still many many temporaries: one for T and one for each matrix product! Also here the matrix T is triangular right ? So I'd suggest the following: Derived::PlainObject tmp(Vtri.transpose() * A); tmp.noalias() = T.transpose().triangularView<Upper>() * tmp; A.noalias() -= Vtri * tmp; There are still 2 temporaries, one for T and for tmp, but that's still better than 4. I'm also not 100% sure about the second line, because in theory a triangular matrix times a matrix can be performed in-place but I have to check the implementation. If you have troubles just remove the .noalias(). But this still looks relatively expensive in term of temporaries, there must be something smarter. gael On Tue, Jun 15, 2010 at 4:08 PM, <vincent.lejeune@xxxxxxxxxx> wrote: > Here is my source code, the part of interest is at the top > > (HHUtilities<..>). > > > > The incrimined lines is the 47/48 : with line 47, everything works, but it > > does a copy I would like to avoid. > > With line 48, I do no copy, but it crashs. > > > > "Derived" is Matrix<float,Dynamic,Dynamic> and OtherDerived is > > Block<MatrixType> where MatrixType is the previous > > Matrix<float,Dynamic,Dynamic> one. > > > > On Tue, 15 Jun 2010 16:02:04 +0200, Gael Guennebaud > > <gael.guennebaud@xxxxxxxxx> wrote: > >> Hi, > >> > >> could you resend your piece of code in a text file please, because > >> inlined in the email, it is really hard to read ! > >> > >> Regarding the previous issues, what was it ? > >> > >> > >> cheers, > >> gael > >> > >> > >> On Tue, Jun 15, 2010 at 3:42 PM, <vincent.lejeune@xxxxxxxxxx> wrote: > >>> > >>> Hi, > >>> > >>> > >>> > >>> I fixed the issue I had that prevent me from getting a working > >>> > >>> implementation of QR Householder decomposition, and apply Householder > > on > >>> > >>> the left/right function. > >>> > >>> > >>> > >>> However, Im' currently facing another issue : my implementation is > > slower > >>> > >>> than the unblocked version of theses functions, because I'm relying on > >>> the > >>> > >>> "triangularView" template. > >>> > >>> As I understand, matrix.triangularView() does not copy the matrix > > object, > >>> > >>> and provides some comfort. However, the following > >>> > >>> applyBlockHouseholderOnTheLeft function does trigger a segfault : > >>> > >>> > >>> > >>> template<typename Derived, typename OtherDerived = Matrix< > >>> > >>> typename Derived::Scalar, Derived::ColsAtCompileTime, Dynamic> > > >>> > >>> class HHutilities > >>> > >>> { > >>> > >>> public: > >>> > >>> typedef Matrix<typename Derived::Scalar, > >>> Derived::ColsAtCompileTime, > >>> > >>> Derived::ColsAtCompileTime> TriangularFactorType; > >>> > >>> typedef typename ei_plain_diag_type<Derived>::type CoeffVect; > >>> > >>> > >>> > >>> TriangularFactorType > >>> > >>> static > >>> > >>> BlockHouseholderTriangularFactor(const Derived& V, > >>> > >>> const CoeffVect& coeffs,const Derived& VV) > >>> > >>> { > >>> > >>> const int k = V.cols(); > >>> > >>> ei_assert(V.cols()==k && V.rows()>=k && "Vector storing Matrix > >>> > >>> must have same columns than coeffs size and more rows than columns"); > >>> > >>> TriangularFactorType ret = TriangularFactorType::Zero(k, k); > >>> > >>> const int j = V.rows(); > >>> > >>> > >>> > >>> for (int i = 0; i < k; i++) > >>> > >>> { > >>> > >>> ret.block(0, i, i, 1) = - coeffs(i) > >>> > >>> * V.block(i, 0, j - i, i).transpose() * V.block(i, i, j > > - > >>> > >>> i, 1); > >>> > >>> ret.block(0, i, i, 1) = ret.block(0, 0, i, i) * ret.block(0, > >>> > >>> i, i, > >>> > >>> 1); > >>> > >>> ret(i, i) = coeffs(i); > >>> > >>> } > >>> > >>> > >>> > >>> return ret; > >>> > >>> } > >>> > >>> > >>> > >>> static > >>> > >>> void > >>> > >>> applyBlockHouseholderOnTheLeft(OtherDerived& A, const Derived& V, > >>> > >>> const CoeffVect& coeffs) > >>> > >>> { > >>> > >>> const TriangularView<Derived,UnitLower>& Vtmp = V.template > >>> > >>> triangularView<UnitLower> (); > >>> > >>> > >>> > >>> const TriangularFactorType& T = > >>> > >>> BlockHouseholderTriangularFactor(Vtmp, > >>> > >>> coeffs,V); > >>> > >>> A -= Vtmp * T.transpose() * Vtmp.transpose() * A; > >>> > >>> } > >>> > >>> } > >>> > >>> > >>> > >>> With OtherDerived=Block<MatrixType>, and > >>> > >>> MatrixType=Matrix<float,Dynamic,Dynamic>. > >>> > >>> I'm forced to make a copy of the TriangularView : > >>> > >>> > >>> > >>> const Derived& Vtmp(V.template triangularView<UnitLower> ()); > >>> > >>> > >>> > >>> (instead of > >>> > >>> const TriangularView<Derived,UnitLower>& Vtmp = V.template > >>> > >>> triangularView<UnitLower> () > >>> > >>> ) > >>> > >>> > >>> > >>> to get my code to work. However, this is a bad workaround, as it makes > > an > >>> > >>> additional copy. > >>> > >>> The overhead is not that big (0.1 s on an old P4 3Ghz for a 256x256 > >>> > >>> matrix), but it should not be there, as theoricaly, blocked computation > >>> are > >>> > >>> better on big matrix...So I hope to improve things if I can get rid of > >>> this > >>> > >>> avoidable matrix copy. > >>> > >>> > >>> > >>> Thanks, > >>> > >>> Vincent > >>> > >>> > >>>

**Follow-Ups**:**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>***From:*vincent.lejeune

**References**:**[eigen] Matrix product between some MatrixType and triangularView<UnitLower>***From:*vincent.lejeune

**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>***From:*Gael Guennebaud

**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>***From:*vincent.lejeune

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>** - Next by Date:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>** - Previous by thread:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>** - Next by thread:
**Re: [eigen] Matrix product between some MatrixType and triangularView<UnitLower>**

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