Re: [eigen] patch for tridiagonal matrices in eigen |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] patch for tridiagonal matrices in eigen*From*: "Gael Guennebaud" <gael.guennebaud@xxxxxxxxx>*Date*: Fri, 2 Jan 2009 00:13:34 +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=LgzpCmgFeDySP7XLkdYfQEms8sn7xI3cssNnAtmDKLg=; b=D6AAoC1Q2oG5f9fsgR0S4WlBscvwWdWRrej3Ko0G1apE1oaaBC3evF5a8CzoEggnKR pidLAzuTX+mP5J/QHqRA+0Tu1mPB4nE05ATBRik7H+ZW6p0vp9oXkvTa2DUP6IyEpQQd +uwS5YbDPskONLp4c6lo3nON94ZxmGLlSAGw8=*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=RlyVJ8uXRweSEOpjhtZHSa500z0R/yoIlldeleQn4/bfWkVvw4lvdcQAemf1uioLFj hSZNTM+GuqIDHYrY6ETj7a0hdScKP8mkUhRdttT2o62m9skVBdEnKj7XyPSO946AMWbe ZwWb7sc5xvDZMbjyoiSFJ258h75SUsrFCHpOM=

Hi Benoit, I mostly agree with your proposal. Two things: 1) I would suggest the following declaration: template<typename Scalar, int Size, bool AsXpr = false> DiagonalMatrix; the main reason: this is easier to use as a "diagonal matrix with storage" class, and this does not prevent to use it as an expression via the additional template parameter. Moreover, such an interface is more suited for a TridiagonalMatrix class or for any other matrices with skyline storage. hm... this last remark let though that all these special matrices could be handled by a unique and generic SkylineMatrix class: template<typename Scalar, int Size, int Bands, bool AsXpr> SkylineMatrix; In practice, DiagonalMatrix could inherit SkylineMatrix<Scalar,Size,1>.... 2) In order to guarantee a high consistency, I'm not 100% sure that DiagonalMatrix should not inherit MatrixBase. It is very important than any function which takes as input an Eigen matrix expression, can also takes as input a DiagonalMatrix. Also, inheriting MatrixBase could avoid a couple of code duplication. For instance, trace() is simply .diagonal().sum(), so once we have overloaded both diagonal() (to return a vector) no need to reimplement neither sum() nor trace... However this would require to make sure that all calls to a MatrixBase<> object are prefixed by ".derived()" that is currently not the case. A robust and general solution to that problem is to go back to pairs of publicfunctionname() _implementationfunctionname() where the former would be implemented only in MatrixBase as: publicfunctionname() { return derived()._implementationfunctionname(); } and MatrixBase() would provide a default implementation of _implementationfunctionname(); Perhaps, one or two macro would be good here. As usual, the drawback is compilation time. Cheers, Gael. On Tue, Dec 30, 2008 at 12:23 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote: > Hi Mauro and Gael, > > Finally I started actually thinking about this discussion (until now > my attitude was "noooo resist entering into yet another technical > discussion you don't have time") > > First let's repeat one comment about the DiagonalMatrix class: Its > coeff() method does an if(). So it's slow (except in the fixed-size > cases where the if() can go away when loops are unrolled, but that's > not the general case). In 99% of practical use cases, one doesn't want > and one doesn't need to pay that price. For example, the diagonal > product is rewritten in such a way that this coeff() method is > bypassed. > > That makes me think that the DiagonalMatrix expression was a bad idea > (my bad, it dates back to early eigen2 days) in the first place. Right > now I think it's completely useless, read on to see what I propose to > have instead. > > First let's look at what happened since this DiagonalMatrix class was written: > 1) Gael wrote the Diagonal product code that actually completely > bypasses the DiagonalMatrix class (i.e. the coeff() method). In Gael's > DiagonalProduct, the DiagonalMatrix is just used as a trivial wrapper > around the diagonal vector (which is directly used). > 2) Mauro is now puzzled as to how to implement TridiagonalMatrix > because the wrong design of DiagonalMatrix is even worse when applied > to tridiagonal matrices. > > Before continuing, let's state a rule to determine when an expression > class (here, DiagonalMatrix) is a bad idea (of course it's a > sufficient condition, not a necessary condition): > > An expression class MyExpression is a bad idea whenever the > MyExpression::coeff() method is slow and in practice it can be > bypassed when implementing actual algorithms. > > For exemple, as we said above, in the case of DiagonalMatrix, the slow > coeff() method is bypassed when implementing the diagonal matrix > product. > > What do I propose to do then? > > If MyExpression represents a special class of matrices for which the > storage can be optimized and certain algorithms can be optimized, then > instead of an expression class MyExpression, I think that what we want > is a special matrix class MySpecialMatrix that does NOT inherit > MatrixBase, so it's not at all an expression. MySpecialMatrix should > have only the storage that it needs, and should be where we implement > the special optimized algorithms. > > For example the DiagonalMatrix class, rethought as a "special matrix", > would look like: > > template<typename DiagonalVector> > class DiagonalMatrix > { > public: > template<typename OtherDerived> > DiagonalProductReturnType<blabla> operator*(const > MatrixBase<OtherDerived>& other); > > protected: > DiagonalVector m_diagonal; > }; > > Here it remains to be seen what DiagonalProductReturnType<blabla> is. > I think that we can make it a plain expression (i.e. inheriting > MatrixBase) because there the coeff() method can be implemented > without an if(). Not sure though, how that plays with vectorization > (depending on the matrix storage order). For TridiagonalProduct > though, it's not possible to make it return an expression (again > because the coeff() would do an if()). Probably we want some meta > logic here, like we have in Product.h, to distinguish cases where we > want an expression and cases where we don't. > > You probably are going to ask: am I talking about a solution for the > diagonal-matrix-expression or for the diagonal-matrix-storage? The > answer is that the DiagonalMatrix above works for both! If you take > DiagonalMatrix<VectorXd> you have a diagonal-matrix-storage. On the > other hand, the diagonal-matrix-expression returned by > MatrixBase<Derived>::asDiagonal() would be a DiagonalMatrix<const > Derived&>. The user is not expected to construct that himself, he just > calls asDiagonal(). > > Another aspect of my proposal is that Product.h does no longer have to > dispatch all the special products, they are now implemented in the > special matrix classes. So it's more extensible and potentially > lighter on the compiler. > > What do you think? > > Cheers, > Benoit > > --- > > ---

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] two technical points: WithAlignedOperatorNew and std::complex casting** - Next by Date:
**[eigen] a few more points...** - Previous by thread:
**Re: [eigen] two technical points: WithAlignedOperatorNew and std::complex casting** - Next by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen**

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