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*: Fri, 19 Dec 2008 12:47:10 -0500*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=a/nRnrKmUv586kJo9yK76I1QLrkinNq6D+7jYAb8PYM=; b=eBS4TQ2BUK00/ltf8x+BG5ZN0mRV+Sgh26Gyw04zw95udJl12cGhh1/Yg6XUtjQlV3 lLGRup0JAEFTKR0wsbyl+HvdGvl4sH/pdSZnYOvUazuiW6b9bken8F0uVR9pL40EFpQ2 6XrHUqAQvb3RaUOS52ycMR1JXmF+wvUqO7lGo=*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=TqOXLIOHY2Zxxp7H4hLhsr817zOzb/izpBs2yYVGyNAppqICpJc3GjZ8A5FBDum70d ip1gH3+vwiXY6TIUW5JkT0Zv6rU+Ipq/zSwla1BxYECOqr1j13PPjc4OKbd4brgy5PX4 Yk44CxBFShV4QkCzS5LVNDh949gdqnlwnTETc=

Hi Mauro, Thanks for the update. I'll let Gael keep commenting on the technical points as he does that much better than me, but here's my take on the release-management part: 1) The Upper ---> UpperTriangular move should happen as soon as possible in order to let people some time to adjust to this API change. So if you can send me a patch with just that, i'd like to include it in beta3 this weekend. Yes it affects a lot of places indeed. 2) About whether to add this to Core or add a new "experimental" area, i'd suggest instead to just wait until we branch 2.0, which should happen in januay, in less than 1 month, at the same time that we release eigen-2.0-rc1. Then trunk will be wide open for invasive commits like this. It's true that currently we're in heavy pre-2.0 bugfixing mode and it's not the best time for such a change until we branch. Until now a lot of people use trunk to build trunk/KDE but that's going to change, see the "kdesupport rearrangement" thread on release-team@kde. So i have good hope that after the branching we can freely break things in trunk. Cheers, Benoit 2008/12/19 Mauro Iazzi <mauro.iazzi@xxxxxxxxx>: > Ok, I'm sorry for the long time it took and thanks for the comments. I > address them below. > > Following what Benoit said I changed Upper to UpperTriangular and now > UpperTriangular==UpperTriangularBit > Unfortunately this affects a lot of places so I think it deserves a > patch on its own. I will check again before submitting. In any case I > agree with Gael about being cautious over the patch because a lot of > things are still unclear to me and I may have missed the proper way of > doing something. > > I do not understand properly the meaning of these flags. For example > part needs assumes some to be exclusive but not others (i.e. > UpperTriangular is not UpperHessenberg, but StrictlyUpper is > UpperTridiagonal...). If part() checked using == instead of & I think > it would be safe to use the Flags to mark the type of the matrix (it > already does in some points). > > This would automatically allow to skip a QR step if the matrix is > already StrictlyUpper just by checking for UpperTriangular... Or > tridiagonalizing a diagonal would be O(1) just checking for a > Tridiagonal (which was why I assumed it worked like that). What do you > think? > > 2008/12/15 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>: >> >> Hi, >> >> wow that sounds really great. I have not much to say. >> >> About the shortcuts, yes Upper must remain unchanged (Upper == >> UpperTriangularBit), and I'm also OK to rename it to UpperTriangular. >> >> About the vectorization of DiagonalProduct, actually some cases was already >> vectorized for a long time. Also note that it does not have the >> EvaluateBeforeAssign bit, so it can already works in-place. What is missing >> is the vectorization of the other cases with all the mess we discussed last >> week... Now we have a TridiagonalProduct as well, I also think that's better >> to go for a C-style function with in-place evaluation. >> >> Marco: in ei_traits<Product<TridiagonalProduct> >::Flags you should remove >> the VectorizableBit flag, we will vectorize it later using an auxiliary >> C-style function. > > you are right: done > >> Also, I think we need two different implementations of TridiagonalProduct: >> one when we take the tridiagonal part of a complete matrix (or when we mark >> a complete matrix as tridiagonal), and another one when the tridiagonal >> matrix relies on a sparse storage (set of 3 vectors). Currently it relies on >> dynamic branching in TridiagonalProduct::coeff(). That's ok for fixed-size >> because our meta-unrollers will allow the compiler to remove them, but >> that's pretty bad for large sizes. Actually, an easy and generic way to >> implement it would be to use a sum of 3 diagonal products, but that would >> not be cache friendly, so... >> >> I also think the TridiagonalMatrix expression should overload the diagonal() >> function (to return a vector) and provide the respective sub-diagonal access >> functions. > > yes, I will do that. > >> We probably also need to add a specialization of assign to optimize simple >> operations on tridiagonal matrices (sum, etc.), more like we did for >> symmetric matrices. > > It would be very good. I will look into it. Making a new Matrix type > and specializing this operations should be easy to do. > >> Then there is the story of tridiagonal symmetric matrices.... >> >> We should also come up with a mechanism to define TriagonalMatrix with their >> own storage. For instance with DiagonalMatrix, one can currently do: >> >> typedef DiagonalMatrix<NestByValue<VectorXf> > DiagonalMatrixXf; >> >> to almost get a "real" diagonal matrix of float with dynamic size. I say >> "almost", because we cannot set its size, resize it, etc... Well one can add >> take a vector, and add a ".asDiagonal()" all the time, but, 1) that's not >> super convenient, and 2) this would be even less convenient for >> TridiagonalMatrix, so we really need to come up with something here. > > I do not follow. What's wrong with DiagonalMatrix<VectorXf>? > >> To conclude, I'm obviously 100% OK that you commit your stuff whenever you >> can, but I'm not in favor to add such new beta features in Core. So I would >> suggest to put the support for tridiagonal matrix in the Sparse module or >> into a new module that we could call Experimental to give us the time we >> need to refine it a bit. What do you think ? > > given how trivial it was to create this patch I think it is not worth > to commit it yet. Anyone in a strong need for tridiagonals can do it > by himself or take this from the ml. in the tree it means just more > code to mantain while at present I'm really not that sure I will have > time. > >> About the QR module, I think that MatrixType is already square, isn't it ? >> so I don't really follow you... About the "if", don't bother too much since >> this need to be redesigned anyway (see the other thread), but your use case >> (the matrix is already tridiagonal) will help to make the right choices... > > Yes, but is particular cases it is not assignable: > SelfAdjointEigenSolver< TridiagonalMatrix<VectorXd, VectorXd, VectorXd> > > fails because it defines m_eivec as a triangular matrix. But if you > already have a tridiagonal matrix you do not want to cast it into a > square matrix before. > > not counting that then the solver cannot understand that it can skip > tridiagonalization unless is checks (n-1)(n-2) coefficients, while it > could get that info from the compiler... > > The SquareMatrixType is the corresponding square type with assignable > coefficients. so I think it is the correct type for m_eivec. Same idea > could apply elsewhere. > > I hope that makes sense. > > cheers, > > mauro > >> >> cheers, >> Gael. >> >> >> >> On Mon, Dec 15, 2008 at 2:59 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> >> wrote: >>> >>> >>> -const unsigned int Upper = UpperTriangularBit; >>> >>> -const unsigned int StrictlyUpper = UpperTriangularBit | ZeroDiagBit; >>> >>> -const unsigned int Lower = LowerTriangularBit; >>> >>> -const unsigned int StrictlyLower = LowerTriangularBit | ZeroDiagBit; >>> >>> +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit; >>> >>> +const unsigned int StrictlyUpper = Upper | ZeroDiagBit; >>> >>> +const unsigned int Lower = LowerTriangularBit | LowerHessenbergBit; >>> >>> +const unsigned int StrictlyLower = Lower | ZeroDiagBit; >>> >>> >>> >>> here I don't really follow. These are parameters for part(). So unless >>> >>> you want to add Hessenberg support in part (which would be a good >>> >>> idea) no need to change that. If you want to add Hessenberg support, >>> >>> you need to make the corresponding changes in Part.h. >>> >> >>> >> I saw that the product helper checked for diagonals using the flags, >>> >> so I did the same for tridiagonals. >>> > >>> > OK right, sorry. These are always convenient shortcuts. What needs to >>> > be updated then is the comment saying that these are the parameters >>> > for part(), actually they're not used only for part(). >>> >>> Wait, no, I answered too fast. We really can't do such a thing as >>> >>> +const unsigned int Upper = UpperTriangularBit | UpperHessenbergBit; >>> >>> for example how would we then in part() distinguish between upper >>> triangular and upper hessenberg... >>> >>> the real problem is that in the first place we allowed ourselves this >>> silly shortcut Upper for upper triangular. Now that there is upper >>> hessenberg too, this is getting ambiguous. I think we just rename >>> Upper to UpperTriangular and then we add UpperHessenberg alongside. >>> >>> Yes this is an incompatible API change, there's a reason why we >>> haven't released yet. KDE / KOffice is not affected, as far as I can >>> remember. >>> >>> 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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] patch for tridiagonal matrices in eigen** - 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/ |