Search code examples
c++matlabeigenmex

How can I build an Eigen::Map to a matlab::data::Array object that I intend to populate with data


This question is a follow up to a previous question about building an Eigen::Map to an input matlab::data::Array object in the Matlab Data API for C++. Now I want to know how I can implement something similar but to an output array which I will populate with data as I loop through it.

In the previous question, I implemented the linked solution from Matlab's forums. Here, I tried using the following to get a non-const pointer:

template <typename T>
T* getOutDataPtr(matlab::data::Array arr) {
  matlab::data::TypedArray<T> arr_t = arr;
  matlab::data::TypedIterator<T> it(arr_t.begin());
  return it.operator->();
}

Which is used like so:

// Initialize input parameters
long int numEl = inputs[0][0];
long int na = inputs[0][1];
long int nPoints = inputs[0][2];
        
// Preallocate output and define map to output array
outputs[0] = factory.createArray<double>
  ({static_cast<size_t>(nPoints),static_cast<size_t>(numEl*na)});
        
auto ptrRecon = getOutDataPtr<double>(outputs[0]);
Eigen::Map<Eigen::MatrixXd> Recon(ptrRecon,nPoints,numEl*na);

for(int i = 0; i < nPoints; i++) {
 for(int n = 0; n < na; n++) {
  for (int j = 0; j < numEl; j++) {
   Recon(i,n*numEl + j) = 5.0;
  }
 }    
}

Even with a simple assignment like that however, I encounter runtime errors after compiling and Matlab crashes. What exactly is going wrong?

Additional Trials I also tried the following for the inner loop which does seem to work:

for(int i = 0; i < nPoints; i++) {
    for(int n = 0; n < na; n++) {
        for (int j = 0; j < numEl; j++) {
            outputs[0][i][n*numEl + j] = 5.0;

        }
    }    
}

I suspect this means there's some issue with the way I am getting the pointer to the output matlab data array?


Solution

  • This is just a guess because I'm going off the rather mediocre Matlab API documentation and don't have a Matlab license to test it with. So please take this with a grain of salt.

    One difference that I see is that you create a mutable iterator TypedIterator<T> instead of the TypedIterator<const T> used in the Matlab forum. And you do this in a function that takes matlab::data::Array by value. This creates a shared data copy. As I understand it, these are copy-on-write.

    My hypothesis is that either constructing that iterator or using it in a manner that could modify that shared array "un-shares" it, effectively creating a copy. However, that copy would be limited to the life-time of the argument, meaning the function. After the function returns, that becomes a dangling pointer.

    Try changing the function signature to T* getOutDataPtr(matlab::data::Array& arr). Or inline the code into your main function.

    Edit

    Actually, that change wouldn't be sufficient because you construct a TypedArray<T> from the Array. That also counts as a shared copy. Try this instead using getWritableElements:

    template <typename T>
    T* getOutDataPtr(matlab::data::Array& arr) {
      auto range = matlab::data::getWritableElements<T>(arr);
      return range.begin().operator->();
    }