[eigen] discussion on special matrices |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: [eigen] discussion on special matrices*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Mon, 30 Mar 2009 16:56:52 +0200*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:date:message-id:subject :from:to:content-type:content-transfer-encoding; bh=lEAh6DOfeYBm7TFEOEIdkj7daQNRjdN8EbSm+z8aZX8=; b=DhKqanIZNrFpwuee/2LoKg1vRm/35Pbf218/aDyrk3nSCsnzA7pETHIs7JNh9kjmfB azkOINz1Jtx4pEnWTIrPE/LVgZG6o+LV5PooQy5VudaavMmGWcsi7SIJvkvJ3dohTy3a 96LNWlxwb2tBZJLNiKyqUNAO6aYh8jgISHWJ4=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type :content-transfer-encoding; b=UcTOZD6Uk1EY7NQ7SFwt2JKW3LoTiIZAai+PWSg2rNvFh9pWjd3EEkefp5Jg0AEvQm IqdmSfthXd7Xw++m5itLouPzMZBDBmzGxOqgzvK+OJFVZbsSK2LUBZK4sYRvH1mKvnxt rTp6M5d3NRFe51Wisjk5EGu/dD/Z1gCBuqpdU=

Hi, After a long break away from coding I am trying to get back to the important discussion on special matrices. I know that's long overdue and hope Mauro didn't lose all hope... So Gael, I studied your page http://eigen.tuxfamily.org/index.php?title=SpecialMatrix and I think I agree with most of your ideas but sometimes I can't be sure because it's hard to tell what you meant. For example you say "remove Part" but then "move solveTriangular* to Part::solve* " so I guess that you meant TriangularView::solve*. Anyway, let me make a proposal following closely yours with a few differences, to see if we agree. I use email for now but once we agree on a solid base, we come back to the wiki. *** CHANGES TO THE PART CLASS *** A. replace if() by ei_assert() in all coefficient accessors in Part. So, to schematize, Part<UpperTriangular>::coeff(int i, int j) { return i<=j ? m_matrix(i,j) : 0; } becomes Part<UpperTriangular>::coeff(int i, int j) { ei_assert(i<=j); m_matrix(i,j) : 0; } In other words, Part becomes completely trivial (useless) from the point of view of coefficient access. The advantage is that it does not have a cost. Notice that we already have this behavior in coeffRef(). B. Part no longer inherits MatrixBase C. We keep the operator= Part::operator=(const MatrixBase&) that allows to do " m.part<UpperTriangular>() = other_xpr;" Yes we'll discuss later the ugliness of the triangular/rectangular mismatch in the API here... I don't see an obvious solution. Maybe we can name it "copyFrom" instead of "operator=" ? D. We add a few methods like setZero() as that is widely used, trivial wrappers around C. E. We add the triangular solving here as Part::solve(). F. If we implement compact triangular storage as a separate class CompactTriangular, then it seems indeed to make sense to move some methods (like solve) to a TriangularBase, however it's nontrivial, because for high performance, the coefficient accesses need to be done differently in a CompactTriangular matrix. Indeed, since every coeff index by row and column takes 7 arithmetic operations, we don't want to implement the solving with this kind of access! Without vectorization, i'd be confident that we can factor coeff accesses using pointers to rows/columns depending on storage orders. So that's the kind of coefficient access that we should implement in a TriangularBase. But to say everything, between the difficulty of implementing efficient triangular solving and the fact that each access requires 7 ops, I'm not convinced about CompactTriangular. Another thing: if we make Part inherit TriangularBase, then we are restricting Part to triangular parts, so that means more classes, and some meta kungfu in MatrixBase::part() to determine the return type. ***** DIAGONAL MATRICES ***** I agree 100% with your proposal. The removal of the MatrixBase inheritance is completely aligned with my chakras of simplifying the whole scheme and making all these special matrix classes just simple, dumb, specialized and efficient. **** OTHERS **** I assume that for the others too, you meant no MatrixBase inheritance (that was my opinion when we first discussed it so I may be wishful here, correct me if I am). For Band matrices, we really need non-square ones. The most important example is bidiagonalization, that is useful for the SVD. It should work for all rectangular matrices (no reason at all to restrict to square). On the other hand, it's also very important to have the special case of self-adjoint band matrices. For tridiagonal matrices, in most cases they're selfadjoint (used for selfadjoint eigensolving) and then we want to optimize the storage. Then of course they're square. ***** SECTION "How to glue everything together" ************ About the Common base class, I don't understand. You suggest to move operator+ there ; first of all I was not convinced that such an operator would be useful at all, I was more seeing these special matrix classes as holding intermediate data to implement and/or call solving/decomposition algorithms, and I don't see the use case for performing a full addition of band matrices....? Then, assuming we really want operator+, it would be implemented in a completely different way between different special matrix classes, so could not share a common implementation in a base class. About providing specializations of assignment, i completely agree. For the same reason as above, I wouldn't bother at all about CwiseBinaryOp. Cheers, Benoit

**Follow-Ups**:**[eigen] Re: discussion on special matrices***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**[eigen] sse asin implementation** - Next by Date:
**[eigen] Re: discussion on special matrices** - Previous by thread:
**Re: [eigen] sse asin implementation** - Next by thread:
**[eigen] Re: discussion on special matrices**

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