Eigen: Return a reference to a block of a matrix w

What I am asking is a generalization of this question. Specifically, I would like to make a C++ Eigen wrapper around a legacy 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 themselves column 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 make extracting entries of this matrix as easy as possible. In C++14, that looks like this 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 that 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's 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;
}

Specifically, 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 issue is that the compile-time checks on the dimensions don't seem to be working. In the data_getter.h, you may notice that I put a bunch of static_asserts on the dimensions. That may seem like a bit of overkill, but I wanted to ensure that the expressions really are performing compile-time operations so that we can get checks on the dimensions. If they were dynamic expressions, then the sizes would all be -1.

However, despite the fact that all of the static_asserts pass, there doesn't seem to be any compile-time checking on the references. For instance, if we change the following line in the test program

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

into

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

The code compiles, but yields a runtime crash. This seems to suggest that the Ref is discarding dimensions. So how should I be defining these variables?

One idea that may come to mind is to define these return values as auto as well. However, this is explictly discouraged by the Eigen docs because if we end up using the output in a loop, it may lead to the expression being evaluated over and over again. That is the reason why I am using Refs. Also, it just seems like a good idea to explicitly state the size since we know it at compile-time...

So is this a bug in Ref? and what type should I be using for the variables that are being spit out by all of my accessor methods?

In case you miss it in the comments... @ggael says that Eigen Refs do not check dimensions at compile-time.