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*: "Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>*Date*: Tue, 30 Dec 2008 12:23:00 +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=ORDY5zkNqgLgTf2oBUkFYyJCStwR5a/Laj2dT5sCBg4=; b=Mzw2RlQW711Dpnv+09797l6992OwyrqBwCUegeSy7kXh+N6rb+KA6okDZE7ceY6bh8 jvL9KwTuYfuUUBnbVgPcyjBRWjW4w6HoiRUjSDaqBYcOKQTy7IPtOMsS6VMWaZw1szBI BOnjqjj2j5GflBSEIphD+RnEhHRvWUVkCljLQ=*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=iuEmGLrL7GHsdeIo+beMi19ze9O6Za9r3yCNJMo1a50czLVWaqSrszIhifR+v2FtX0 JAb/bF6E3CtFvUe2ReEuuBdJzNawUMtT6krkWcKlmR5HzM2l1YF+2JKWHmTCu8UyL24U v90sYpZO2l8BSBcLiEClaMgLlSbmIu6ougXSk=

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

**References**:**[eigen] patch for tridiagonal matrices in eigen***From:*Mauro Iazzi

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

**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Mauro Iazzi

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

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

**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Gael Guennebaud

**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Mauro Iazzi

**Re: [eigen] patch for tridiagonal matrices in eigen***From:*Gael Guennebaud

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Next by Date:
**Re: [eigen] memory savings in Tridiagonalization::decomposeInPlace** - Previous by thread:
**Re: [eigen] patch for tridiagonal matrices in eigen** - Next by thread:
**[eigen] early beta3-preview, MSVC/SSE2 testers wanted**

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