[eigen] Feature idea: coeffRef() ---> setCoeff()
• 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

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

(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>
//...
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

```

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