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*: "Mauro Iazzi" <mauro.iazzi@xxxxxxxxx>*Date*: Thu, 8 Jan 2009 18:59:18 +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=4alLzNlSr3+vR/Y6JJFvU+GUm/9bh/tqMXzeNCzOhQs=; b=fvm+wLF6MnK8Vcydd8yIkzID+e0zVJnEegtJrU0S0WkDqU8KB0tzbyaFGbnWWv0vxZ GG8lr7ISN8iZ+6arRyse4aJAp0iBf0ADGohK5Me0yW8z5vvGHvP1qUj68N5/zpCC/PgF sVxHJkLdp/EPfAnzGouvL3/KBezdSGO7YAVd8=*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=HLaKoA8ZrKULZi+0P5ed+D7y6i8KCpGDpwpyor7z1nJ6d00YuSD341l699ThVieRil lKDmyAKBDoeVAFA1iPfHKnbtVC+ObYRDzbbGtMm/74vmGG4B1KDep1AS0IOu+Wa5ajRs YnW2AXQdyenTq1Ko5W7Op8+mfKOR8r/nwIgkA=

Hi, sorry for replying so late. I agree with most of Benoit' and Gael's ideas, but I have some comments relative to my use case. Also, I'm starting to understand this issue better, thanks to gael's explanation, but I may have missed something. Tell me if this is the case. 2008/12/30 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>: > > 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. > It's a perfect criterion. However, when something else than an expression is used, memory may be wasted. This is quite a concern with me: a single matrix storage may be the few hundred megs that exceed my memory and start using swap. I have no obvious solutions. > 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; > }; > I think this is a good idea, since someone could implement a custom storage class in this self contained way. However, product by a diagonal matrix can be optimized for both left and right multiplication. Moreover, I would not rely on template specialization to choose the right path (like template<...> ReturnType<> operator+(const MatrixType &a, const DiagonalType &b); ) but this is because I find this things quite unpredictable. As Gael, I think it could still derive from MatrixBase. It seems to me that a lot of things should be reimplemented in the opposite case. > 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. I have to add that the sum shows similar problems. and while the product expression between a diagonal and a square matrix can be written without if()s, the sum expression cannot (unless I am missing some obvious way). So also just for diagonals the problem remains. In general return expressions should follow the general rule for Bad Ideas you stated before. I think it's worth listing a few cases, assuming there is a DiagonalType and TridiagonalType DiagonalType D; TridiagonalType T; MatrixType A; A = D + A; A = T * A; A = D + T; all of those can be optimized if written directly. If one uses an expression there are if()s and if one does not it creates an unnecessary temporary (which can be lethal in my case). Last case has the strongest possible optimization. Things like A = (D + A) * A; make this problem even more so. Seems to me that the problem is one cannot optimize unless he knows where to put the result (e.g. when there is an assigment). One possible solution would be to allow the class to specify if it can do better than the normal algorithm. Or, even better to use a cost system like the one that is already in place for expressions vs. cache. So the return type of (D+A) could override the assignment to the fast version. On multiplication it would create a temporary as normal. (D+T) could return a tridiagonal expression etc... I guess this is just a more complicated version of what Gael suggests. So even slower compilation times, since a lot of paths would need to be specified. I can't really say this is a good solution. But it seems one that saves memory. A less complicated and less user-friendly method would be to just use expressions and allow for specialization of eval() (and if possible of assignment). In this way one can choose what he wants to optimize. A relatively sane default for these expressions' flags can make this solution quite painless in the general case. > 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. I'd be really happy about this. I recently had a sparse matrix defined in a way that makes it very easy to compute the product without even saving the coefficients. If it was possible to define an expression for it cheaply it would be a big gain. I hope this made sense. Cheers, mauro ---

**Follow-Ups**:**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] State of the problem with std::vector** - Next by Date:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Previous by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Next by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen**

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