[eigen] Feature idea: coeffRef() ---> setCoeff() |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: [eigen] Feature idea: coeffRef() ---> setCoeff()*From*: Manoj Rajagopalan <rmanoj@xxxxxxxxx>*Date*: Thu, 29 Jul 2010 19:39:02 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

Hi eigen developers, Grepping for coeffRef() in src/Core/*.h seems to show that all occurrences seem to be returning Scalar& for one of two purposes: (1) for use as lvalue in assignments, (2) in return expressions within member functions so that users may do the same at a higher level. If this is indeed the case, then I think replacing this with the following will allow greater extensibility: setCoeff(row,col,value) in short, we can have expressions like: MatrixXcd A(N,N); A.conjugate().setCoeff(i, j, value); A.selfadjointView<Lower>()(j, j+1, value); // note: write to upper triangle The nice '=' syntax is not present at the coefficient-wise level but, * it is still present at matrix-expression level, * but we gain a few more advantages (see below) In some detail: 1. The Transpose<> class has been created to allow writes to the nested MatrixType even though the structure has been changed. This is because the Transpose class *logically* represents the structure transformation in a *reversible* way so as to allow coherent definitions of coeff() and coeffRef(). 2. This idea can be extended to conjugates. A Conjugate<> wrapper, like Transpose<>, will reversibly encapsulate the coefficient-wise conjugate operation. IMPORTANT: it will enable writes because conjugation is a bijection on real/complex fields but this will depend on the existence of either of two alternatives: (a) a setCoeff() method as opposed to a "coeffRef() =" kind of assignment which reverse-transforms values before writing to memory, OR (b) a wrapper class for Scalar& that is returned in coeffRef (instead of Scalar&), whose operator=() is overloaded to perform the reverse-value-transformation. 3. The benefit here is that we may define a HermitianMatrix or SelfAdjointMatrix class for which writes to either triangle will now be coherent from the user's view: Eigen can internally map the write down to the stored/viewed triangular portion (I am working on CompactSelfAdjointMatrix and related classes that do this - will release soon). (a) We may even define SkewHermitianMatrix where the elements must be conjugated and negated between the two triangles. Not possible easily with existing architecture. (b) Such specialized types capture important physical information about the problem at reduced memory requirements and improve bug-safety. 4. Generally, in any view of an underlying dense matrix type, coherent writes to the coefficients can be performed if the view transforms each coefficient in a 1-1 lossless manner. Structural rearrangements like transpose and permutations allow direct access of Scalar& because only a computation on the indices is involved. But value-transformations like conjugate() will require a setCoeff() method because the value is transformed (reversibly). 5. By allowing both structural and value reversible transformations simultaneously, we can extend Eigen by composing these features. For example, the return type of adjoint() can now be Conjugate<Transpose<MatrixType> > which can now *allow writes* through the setCoeff() call-chain. I imagine this will be useful in some inplace-solve functions. Currently, AdjointReturnType is built similarly, but because conjugate() is a "const" method, no writes are possible. 6. FAR-FETCHED OPTIONAL IDEA: The above idea can be extended to other reversible value-transformations if needed (by the user). For example, the following contrived example becomes possible for free: VectorXd v(N); (5.0*v).setConstant(1.0); // v[i]=0.2 for all i This form of expression might be natural in some contexts where scales are important (eg. performing energy computations in eV instead of J). (a) why the user would do the above is upto the user but the functionality is available for free because ei_scalar_product_op is 1-1 for non-zero scales (b) For many ops like scalar_sum_op and scalar_difference_op, the "inverse" op is easy to derive and in fact can be templated in the traits class. Eg: template<typename Scalar> struct ei_traits<ei_scalar_add_op> { //... typedef ei_scalar_subtract_op inverse_op; // user has to instantiate with negative of add-op's parameter though }; If for each op, an inverse op can be so defined then the architecture will naturally allow synthesis of write-able cwise-transformed views of dense matrix storage types. (c) it would be up to the user to avoid non-sensical parameters like 0 for the product which would render the inverse undefined (d) (5.0*v) is a lossy transformation: bits are lost in finite-precision arithmetic - the losslessness is only in the algebraic definition of the * operator but as long as the user is OK with this, he/she can use these features. 7. I've included the above far-fetched idea as an example of how users can extend features beyond what the developers intended. Template Metaprogramming itself is a result of such serendipity. The above example might be useless but who knows someone might figure out something really useful :-) Thanks, Manoj

**Follow-Ups**:**Re: [eigen] Feature idea: coeffRef() ---> setCoeff()***From:*Benoit Jacob

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] licensing: watching the MPL update process** - Next by Date:
**Re: [eigen] licensing: watching the MPL update process** - Previous by thread:
**Re: [eigen] licensing: watching the MPL update process** - Next by thread:
**Re: [eigen] Feature idea: coeffRef() ---> setCoeff()**

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