[eigen] Eigen Types as Parameters for Functions

[Christian Seiler <christian@xxxxxxxx>]

Hi,

I'm a user and a fan of Eigen. However, there is one issue I frequently
run into: I very often want to write functions that take Eigen types as
parameters. Now, I've read
<http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html>,
and that has certain issues: First of all, these functions have to be
defined as template functions, which means that they have to be defined
in a header file in order to work. This is fine if the function does 
not
do that much, but I'd really prefer to put the function into a separate
compilation unit if the code is a little more complex. (But for this to
work with template functions, I'd have to explicitly instantiate all
possible template variants in the compilation units, and there I'm
definitely going to miss quite a few.) Also, having to play around with
const_cast<> and template logic seems a bit excessive to me just to
pass a simple matrix / vector as an argument, while at the same time
being able to accept Block-type objects or direct expressions etc.
Especially if the code is supposed to be shared with other people that
are not as well-versed in the intricacies of C++ as I am, who may also
need to modify the code later on.

Therefore, I'd like to propose a solution that would make the life for
me (and probably a lot of other people) easier.

Rationale: In the end, the only real information about a matrix stored
in memory that one needs are:
    - pointer to the first element
    - dimension
    - row/col major
    - strides (inner/outer)
    - alignment

I propose adding a new class (name could be decided on later, for the
purpose of this class, I'll use MatrixParameter, although I don't like
the name particularily) that has the following properties:

 - subclass of MatrixBase, i.e. can be used inside the function like a
   Matrix object
 - has private members containing the above-mentioned bits of
   information (pointer to first element, dimension, ...)
 - can be implicitly converetd to from any compatible MatrixBase (i.e.
   foo(mat) and foo(mat.block(...)) can be used directly)

This would allow one to create functions and/or methods that accept
Eigen objects as parameters much more easily - and would even provide
a stable binary ABI for libraries one wants to write, as long as the
numerical type is specified.

Examples of possible usage pattern:

void veryComplicatedAlgorithm(
  Eigen::MatrixParameterXd result,
  Eigen::MatrixParameterXd orig
) {
  // do some very complicated calculation with orig
  // and store it in result
}

Matrix4d orig;
orig <<  1,  2,  3,  4,
         5,  6,  7,  8,
         9, 10, 11, 12,
        13, 14, 15, 16;
MatrixXd result(2,2);
veryComplicatedAlgorithm(result, orig.block(1,1,2,2));

double worksOnlyWithTwoByTwo(Eigen::MatrixParameter2d inp)
{
  // or something...
  return inp(0,0) * inp(1,1) / (inp(0,1) + inp(1,0));
}

Matrix4d foo = Matrix4d::Identity();
MatrixXd foo2 = MatrixXd::Identity(4, 4);
double p = worksOnlyWithTwoByTwo(foo); // compile-time error
double p2 = worksOnlyWithTwoByTwo(foo2); // run-time error
                                         // unless range-checking
                                         // disabled
double p3 = worksOnlyWithTwoByTwo(foo.block(0, 0, 2, 2)); // works
double p4 = worksOnlyWithTwoByTwo(foo2.block(0, 0, 2, 2)); // works

void someWrapperAroundLegacyFortranCode(
  Eigen::MatrixParameterXd mat,
  Eigen::VectorParameterXd vec
) {
  // or something
  xxyyzz_(mat.size(), mat.data(), vec.data());
}

Considerations:

 - The layout of the fields in the class should be fixed once and for
   all, this allows binary compatibility even when using a different
   version of Eigen for compiling both compilation units.
   ("once and for all" should mean something like "only change this
   if the major version of Eigen changes")

 - If the functions are not templated themselves, this fixes the scalar
   data type. I don't see this as a problem though, because that may
   be what people want. Also, this scheme still allows for declaring
   something along the lines of template<typename Scalar> void
   foo(Eigen::MatrixParameter<Scalar, ...> ...); if this is really
   wanted

 - Resizing of output matrices is not easily possible here. However,
   I don't think this is necessarily a problem, since most of the
   time the calling code will already know what the dimension of
   matrices used for output have to be. And if somebody really
   requires the ability to resize matrices, it should always be
   possible to use the current way of passing Eigen matrices.

 - Expressions such as foo(A * B) can be tricky. I propse that we
   define two classes: The second class could be prefixed with "Const"
   (e.g. MatrixConstParameter). There, any expression would be
   evaluated into a temporary before being passed to the function,
   while plain matrices and .block(...) would remain untouched. The
   non-"Const"-variant, on the other hand, would only accept plain
   matrices and .block(...) expressions, since only those can
   guarantee write acccess. (Ok, .transpose() would also work because
   one can just flip the RowMajor/ColMajor bit, but .adjoint() would
   already be problematic.)

 - The documentation should warn about memory management caveats here.
   Basically, since the classes proposed here are just wrappers for
   pointers, if the matrix pointed to falls out of scope or is
   deleted, and somebody stored a reference somewhere, this will
   lead to bugs. Unfortunately, disabling the copy constructor is
   not possible, since it is needed to be able to pass through
   arguments to other functions.

I'd like to here your thoughts on this, and I'm offering to implement
my proposal if you agree with the basic idea.

Thanks,
Christian