All geek questions in one place

Eigen: returns a link to a matrix block with compile-time measurement checks

What I ask is a generalization of this question . In particular, I would like to create a C ++ Eigen wrapper around the old C and Fortran library, which uses a 2D data structure:

[ x[0,0] ... x[0,w-1] ] [ u[0,0] ... u[0,w-1] ] [ ... ] [ x[c-1,0] ... x[c-1,w-1] ] [ u[c-1,0] ... u[c-1,w-1] ]

where each of the entries x[i,j] and u[i,j] are the very columns of vectors of size ( nx1 ) and ( mx1 ) respectively. This leads to some complicated (and error prone) pointer arithmetic, as well as some very unreadable code.

Therefore, I want to write an Eigen class whose sole purpose is to simplify the extraction of records of this matrix as much as possible. In C ++ 14, this looks like data_getter.h :

 #ifndef DATA_GETTER_HEADER #define DATA_GETTER_HEADER #include "Eigen/Dense" template<typename T, int n, int m, int c, int w> class DataGetter { public: /** Return a reference to the data as a matrix */ static auto asMatrix(T *raw_ptr) { auto out = Eigen::Map<Eigen::Matrix<T, (n + m) * c, w>>(raw_ptr); static_assert(decltype(out)::RowsAtCompileTime == (n + m) * c); static_assert(decltype(out)::ColsAtCompileTime == w); return out; } /** Return a reference to the submatrix * [ x[i,0], ..., x[i,w-1]] * [ u[i,0], ..., u[i,w-1]] */ static auto W(T *raw_ptr, int i) { auto out = asMatrix(raw_ptr).template middleRows<n + m>((n + m) * i); static_assert(decltype(out)::RowsAtCompileTime == (n + m)); static_assert(decltype(out)::ColsAtCompileTime == w); return out; } /** Return a reference to the submatrix [ x[i,0], ..., x[i,w-1]] */ static auto X(T *raw_ptr, int i) { auto out = W(raw_ptr, i).template topRows<n>(); static_assert(decltype(out)::RowsAtCompileTime == n); static_assert(decltype(out)::ColsAtCompileTime == w); return out; } /** Return a reference to x[i,j] */ static auto X(T *raw_ptr, int i, int j) { auto out = X(raw_ptr, i).col(j); static_assert(decltype(out)::RowsAtCompileTime == n); static_assert(decltype(out)::ColsAtCompileTime == 1); return out; } /** Return a reference to the submatrix [ u[i,0], ..., u[i,w-1]] */ static auto U(T *raw_ptr, int i) { auto out = W(raw_ptr, i).template bottomRows<m>(); static_assert(decltype(out)::RowsAtCompileTime == m); static_assert(decltype(out)::ColsAtCompileTime == w); return out; } /** Return a reference to u[i,j] */ static auto U(T *raw_ptr, int i, int j) { auto out = U(raw_ptr, i).col(j); static_assert(decltype(out)::RowsAtCompileTime == m); static_assert(decltype(out)::ColsAtCompileTime == 1); return out; } /** Return a reference to the submatrix * [ x[0,i], ..., x[c-1,i]] * [ u[0,i], ..., u[c-1,i]] */ static auto C(T *raw_ptr, int i) { auto out = Eigen::Map<Eigen::Matrix<T, n + m, c>>( asMatrix(raw_ptr).col(i).template topRows<(n + m) * c>().data()); static_assert(decltype(out)::RowsAtCompileTime == (n + m)); static_assert(decltype(out)::ColsAtCompileTime == c); return out; } /** Return a reference to the submatrix [ x[0,i], ..., x[c-1,i]] */ static auto Xc(T *raw_ptr, int i) { auto out = C(raw_ptr, i).template topRows<n>(); static_assert(decltype(out)::RowsAtCompileTime == n); static_assert(decltype(out)::ColsAtCompileTime == c); return out; } /** Return a reference to the submatrix [ u[0,i], ..., u[c-1,i]] */ static auto Uc(T *raw_ptr, int i) { auto out = C(raw_ptr, i).template bottomRows<m>(); static_assert(decltype(out)::RowsAtCompileTime == m); static_assert(decltype(out)::ColsAtCompileTime == c); return out; } }; #endif /* DATA_GETTER_HEADER */

and here is a test program demonstrating how this works:

 #include <iostream> #include <vector> #include "Eigen/Dense" #include "data_getter.h" using namespace std; using namespace Eigen; template<typename T> void printSize(MatrixBase<T> &mat) { cout << T::RowsAtCompileTime << " x " << T::ColsAtCompileTime; } int main() { using T = double; const int n = 2; const int m = 3; const int c = 2; const int w = 5; const int size = w * (c * (n + m)); std::vector<T> vec; for (int i = 0; i < size; ++i) vec.push_back(i); /* Define the interface that we will use a lot */ using Data = DataGetter<T, n, m, c, w>; /* Now let map that pointer to some submatrices */ Ref<Matrix<T, (n + m) * c, w>> allData = Data::asMatrix(vec.data()); Ref<Matrix<T, n, w>> x1 = Data::X(vec.data(), 1); Ref<Matrix<T, n, c>> xc2 = Data::Xc(vec.data(), 2); Ref<Matrix<T, n + m, c>> xuc2 = Data::C(vec.data(), 2); Ref<Matrix<T, n, 1>> x12 = Data::X(vec.data(), 1, 2); cout << "Data::asMatrix( T* ): "; printSize(allData); cout << endl << endl << allData << endl << endl; cout << "Data::X( T*, 1 ) : "; printSize(x1); cout << endl << endl << x1 << endl << endl; cout << "Data::Xc( T*, 2 ) : "; printSize(xc2); cout << endl << endl << xc2 << endl << endl; cout << "Data::C( T*, 2 ) : "; printSize(xuc2); cout << endl << endl << xuc2 << endl << endl; cout << "Data::X( T*, 1, 2 ) : "; printSize(x12); cout << endl << endl << x12 << endl << endl; /* Now changes to x12 should be reflected in the other variables */ x12.setZero(); cout << "-----" << endl << endl << "x12.setZero() " << endl << endl << "-----" << endl; cout << "allData" << endl << endl << allData << endl << endl; cout << "x1" << endl << endl << x1 << endl << endl; cout << "xc2" << endl << endl << xc2 << endl << endl; cout << "xuc2" << endl << endl << xuc2 << endl << endl; cout << "x12" << endl << endl << x12 << endl << endl; return 0; }

In particular, it produces the following output (as expected):

 Data::asMatrix( T* ): 10 x 5 0 10 20 30 40 1 11 21 31 41 2 12 22 32 42 3 13 23 33 43 4 14 24 34 44 5 15 25 35 45 6 16 26 36 46 7 17 27 37 47 8 18 28 38 48 9 19 29 39 49 Data::X( T*, 1 ) : 2 x 5 5 15 25 35 45 6 16 26 36 46 Data::Xc( T*, 2 ) : 2 x 2 20 25 21 26 Data::C( T*, 2 ) : 5 x 2 20 25 21 26 22 27 23 28 24 29 Data::X( T*, 1, 2 ) : 2 x 1 25 26 ----- x12.setZero() ----- allData 0 10 20 30 40 1 11 21 31 41 2 12 22 32 42 3 13 23 33 43 4 14 24 34 44 5 15 0 35 45 6 16 0 36 46 7 17 27 37 47 8 18 28 38 48 9 19 29 39 49 x1 5 15 0 35 45 6 16 0 36 46 xc2 20 0 21 0 xuc2 20 0 21 0 22 27 23 28 24 29 x12 0 0

The problem is that compile-time runtime checks on dimensions do not work. In data_getter.h you can notice that I have imposed a bunch of static_assert on the dimensions. This may seem a little redundant, but I wanted to make sure that the expressions actually do compilation operations so that we can get the size checks. If they were dynamic expressions, then the dimensions would be -1.

However, although the entire static_assert pass passes, link references are not checked. For example, if we change the following line in a test program

 Ref<Matrix<T, (n + m) * c, w>> allData = Data::asMatrix(vec.data());

in

 Ref<Matrix<T, (n + m) * c + 1, w>> allData = Data::asMatrix(vec.data());

The code compiles, but crashes at runtime. This seems to indicate that Ref discarding sizes. So how should I define these variables?

One idea that may come to mind is to define these return values ​​as auto . However, this is clearly discouraged by the Eigen docs , because if we end up using output in a loop, it can cause the expression to evaluate again and again. It is for this reason that I use Ref s. Also, it just seems like I want to clearly indicate the size, since we know it at compile time ...

So is this a bug in Ref? and what type should I use for variables that spit out with all my access methods?

+5
source share
1 answer

If you skip this in the comments ... @ggael says that Eigen Ref does not check measurements at compile time.

+2
source

Source: https://habr.com/ru/post/1275599/

More articles:


All Articles