Re: [eigen] Eigen Types as Parameters for Functions

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2012/1/26 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)

Just note that this is only possible for those matrix expressions
whose coefficients are stored in memory with a simple enough layout
that can be described in terms of strides. That is what we call a
"Direct Access" matrix expression; the test is ExpressionType::Flags &
DirectAccessBit.

With that said, yes, we have been considering doing this, that's bug 58:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=58

If someone wants to do this and needs mentoring, I should find time to
do at least that mentoring.

Cheers,
Benoit


>
> 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
>
>
>