Eigen: Efficient way of referencing ArrayWrapper

Question

I am interfacing some code with raw pointers. So I have extensive use of the map class:

void foo(T* raw_pointer){
    const int rows = ...;
    const int cols = ...;
    Map<Matrix<T, rows, cols>> mat(raw_pointer);

    // DO some stuff with "mat"
}

Now I want to apply some cwise operations in foo, which I accomplish using .array(). The code works, however, it looks very messy due to all of the .array() calls strewn in the function. For instance, for the sake of argument, let's suppose that the function looked like this:

void foo(T* raw_pointer){
    const int rows = ...;
    const int cols = ...;
    Map<Matrix<T, rows, cols>> mat(raw_pointer);

    for (int i = 0 ; i < 1000 ; ++i)
        ... something = i * mat.row(1).array() * sin(mat.row(4).array()) + mat.col(1).array();
}

Part of the problem with this is that it is very unclear what the code is actually doing. It would be much nicer if gave the variables names:

void foo(T* raw_pointer){
    const int rows = ...;
    const int cols = ...;
    Map<Matrix<T, rows, cols>> mat(raw_pointer);

    Matrix<T, 1, cols> thrust = mat.row(1);
    Matrix<T, 1, cols> psi = mat.row(4);
    Matrix<T, 1, cols> bias = mat.row(2);

    for (int i = 0 ; i < 1000 ; ++i)
        ... something = i * thrust.array() * sin(psi.array()) + bias.array();
}

But it would be even nicer if I could get directly get a reference to the ArrayWrappers so that we aren't making any copies. However, the only way I can figure out how to get that to work is by using auto:

void foo(T* raw_pointer){
    const int rows = ...;
    const int cols = ...;
    Map<Matrix<T, rows, cols>> mat(raw_pointer);

    auto thrust = mat.row(1).array();
    auto psi = mat.row(4).array();
    auto bias = mat.row(2).array();

    for (int i = 0 ; i < 1000 ; ++i)
        ... something = i * thrust * sin(psi) + bias;
}

This code works, and upon testing appears to reference the entries in the pointer (as opposed to making copies like in the previous snippet). However, I am concerned about its efficiency since the Eigen documentation explicitly suggests NOT doing this. So could somebody please what the preferred way to define the types for the variables is in such a circumstance?

It seems to me like I should be using a Ref here, but I can't figure out how to get that to work. Specifically, I have tried replacing auto with

Eigen::Ref<Eigen::Array<T, 1, cols>>

and

Eigen::Ref<Eigen::ArrayWrapper<Eigen::Matrix<T, 1, cols>>>

but the compiler doesn't like either of those.

Two things: The `Ref` doesn't make a copy; You can use a `Map — Avi Ginsburg, Mar 12 '18 at 10:01
Hmmm. The only problem is that we have a general interface class for getting an `Map>(mapMatrix.data())` on that object — bremen_matt, Mar 12 '18 at 10:07
Hmm. It seems that that will indeed work. The only issue I am going to have is the loss of compile-time size checking. I think I can live with that though. You should post an answer... — bremen_matt, Mar 12 '18 at 10:12
I have an interface file (https://stackoverflow.com/questions/49008788/eigen-return-a-reference-to-a-block-of-a-matrix-with-compile-time-dimension-che) for accessing the raw data — bremen_matt, Mar 12 '18 at 10:19
That interface removes a lot of the boilerplate that I would otherwise need to have. It also ensures that the raw data is accessed consistently, so it should help reduces bugs. It also does size checking — bremen_matt, Mar 12 '18 at 10:20
Now if set `Map> mapMatrix = ` output of one of those functions, then I think I could `Map> mapArray(mapMatrix.data());` — bremen_matt, Mar 12 '18 at 10:22
Got it. And I see that [ggael](https://stackoverflow.com/questions/49008788/eigen-return-a-reference-to-a-block-of-a-matrix-with-compile-time-dimension-che#comment85385509_49147375) has since fixed the compile time checks. — Avi Ginsburg, Mar 12 '18 at 10:23
However, in the definition of `mapArray`, I would have to redefine the dimensions. So I could potentially introduce a bug there that cannot be checked at compile-time — bremen_matt, Mar 12 '18 at 10:23
Ahh. Actually, this approach didn't work because of the strides in the underlying pointer when using my interface file. But it would work for the original question... — bremen_matt, Mar 12 '18 at 12:40
I'm not following. The third template parameter for `Eigen::Map` is the stride type. The third or fourth parameter of the `Map` c'tor is a runtime modifiable stride. — Avi Ginsburg, Mar 12 '18 at 13:04
I would have to mess around with that. I guess that would be a viable alternative. But I do not have time to look into that now. — bremen_matt, Mar 12 '18 at 13:18
To get a `Ref` with inner stride, write `Eigen::Ref, 0, Eigen::InnerStride<> >` (this is necessary, if you want to refer to a row of a column-major matrix/array). — chtz, Mar 12 '18 at 17:58

score 3 · Answer 1 · answered Mar 12 '18 at 10:20

3

To avoid having to write array() every time you use the Map<Eigen::Matrix... you can use a Map<Eigen::Array... instead/in addition. This will use the default element-wise operators instead of the matrix operators. To use a matrix operator instead, you can use map.matrix() (similar to what you have in your post mat.array()).

answered Mar 12 '18 at 10:20

Avi Ginsburg

10,323
3
29
56

Similarly, you can/need to write `Array thrust = mat.row(1);` etc. – chtz Mar 12 '18 at 17:53

score 2 · Answer 2 · answered Mar 12 '18 at 10:02

auto thrust = [](auto&&mat){return mat.row(1).array();};
auto psi = [](auto&&mat){return mat.row(4).array();};
auto bias = [](auto&&mat){return mat.row(2).array();};

for (int i = 0 ; i < 1000 ; ++i)
    ... something = i * thrust(mat) * sin(psi(mat)) + bias(mat)

has names. And the array wrappers don't persist.

Eigen: Efficient way of referencing ArrayWrapper

2 Answers2