Изменение переменных, которые не должны быть затронуты (переполнение памяти?)

API-интерфейс mex, хотя и доступен на C++, на самом деле является API-интерфейсом C, и его прямое использование делает код похожим на код C. Это не очень хорошее место.

Итак, давайте сначала напишем код, предполагая, что существуют хорошие оболочки массива MEX C++11.

Начнем с параметров функции. Предположим, что DoubleArray и LogicalArray будут представлять массивы заданного типа и иметь конструкторы, которые инициализируют их из правых параметров или создают новые для левых параметров. Массивы являются просто дескрипторами базовых данных Matlab и представляют собой логический указатель. Они будут увеличиваться, чтобы компенсировать их положение по отношению к базовым данным.

Мы не должны использовать массивы фиксированного размера в стиле C, вместо этого следует использовать тип std::array.

struct Params {
   DoubleArray s, t, X, inGridSize;
   size_t NumPts;
   std::array<size_t, 2> GridSize;
   DoubleArray Points, Jacobian;
   LogicalArray Valid;
   std::array<double, 16> M;

   Params(mxArray * plhs[], const mxArray * prhs[]) :
      s(prhs, 0),
      t(prhs, 1),
      X(prhs, 2),
      inGridSize(prhs, 3),
      NumPts(s.size()),
      Points(1, NumPts, plhs, 0),
      Jacobian(NumPts, X.size(), plhs, 1),
      Valid(NumPts, 1, plhs, 2)
   {
      std::copy(inGridSize.begin(), inGridSize.end(), GridSize.begin());
   }
};

Затем эти параметры могут быть переданы различным функциям. Их можно скопировать, если мы хотим изменить любую из позиций массива.

Это довольно непротиворечиво, хотя можно утверждать, что, возможно, inGridSize не нужен, и мы должны инициализировать GridSize из временного DoubleArray.

Дополнительные параметры для расчета отдельной точки можно аналогичным образом обернуть в структуру:

struct PointParams {
   unsigned int Vpts;
   unsigned int Xlen;
   std::array<double, 16> J;
   std::vector<unsigned int> Crds { Xlen };
   std::array<unsigned int, 16> RelCrds;

   PointParams(const Params & p) :
      Vpts(p.GridSize[1] + 3),
      Xlen((p.GridSize[0] + 3)*Vpts)
   {}
};

Наконец, функция Simple_Separate_Parametrics вычисляет результат типа Parametrics:

struct Parametrics {
   double t_crd_dbl, t_frac;
};

Parametrics Simple_Separate_Parametrics(double t_full) {
   Parametrics p;
   p.t_crd_dbl = floor(t_full);
   p.t_frac = t_full - p.t_crd_dbl;
   return p;
}

Calc_Point_... ссылается на Params и PointParams:

void Calc_Point_and_XJacobian_On_SplurfaceMap(Params & p, PointParams & pp) {
   mexPrintf("Begin Inner Function %u, %u\n", p.GridSize[0], p.GridSize[1]);

   auto pms = Simple_Separate_Parametrics(p.t[0]);

   mexPrintf("End Inner Function %u, %u\n", p.GridSize[0], p.GridSize[1]);
}

Points_and_XJacobian_... принимает свои параметры по значению, копируя их так, чтобы указатели массива могли изменяться (увеличиваться) в цикле. Поскольку p изменяется в цикле, на каждой итерации мы создаем новый экземпляр PointParams на основе p.

void Points_and_XJacobian_on_Splurface(Params p) {
   mexPrintf("Begin Outer Function %u, %u\n", p.GridSize[0], p.GridSize[1]);

   for (unsigned int n = 0; n < 1; ++n) {//NumPts
      PointParams pp(p);
      Calc_Point_and_XJacobian_On_SplurfaceMap(p, pp);
      p.Points += 1;
      p.Valid += 1;
      p.s += 1;
      p.t += 1;
   }

   mexPrintf("End Outer Function %u, %u\n", p.GridSize[0], p.GridSize[1]);
}

mexFunction проверяет правильность подсчета аргументов, инициализирует начальный Params и запускает рабочую функцию Points_and_....

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
   // [Points, Jacobian, Valid] = Function( s, t, X, GridSize)

   if (nrhs != 4)
      mexErrMsgTxt("Need four input parameters");
   if (nlhs != 3)
      mexErrMsgTxt("Need three output parameters");
   if (mxGetNumberOfElements(prhs[3]) != 2)
      mexErrMsgTxt("GridSize needs two elements");

   Params p { plhs, prhs };
   mexPrintf("Begin Mex %u, %u\n", p.GridSize[0], p.GridSize[1]);
   Points_and_XJacobian_on_Splurface(p);
   mexPrintf("End Mex %u, %u\n", p.GridSize[0], p.GridSize[1]);
}

Вы заметите, что ни в одном из этих кодов не происходит явного выделения памяти.

Итак, как нам реализовать волшебный соус? Через общий класс MxArray, параметризованный специализацией Traits, зависящей от типа значения. Конечно, его можно повторно использовать в любом файле MEX.

Тип массива действует как массив C, который знает свои границы. Вы можете легко добавить проверку границ в реализации operator[]. Вы также можете добавить смещение к массиву, сместив позицию элемента на 0 индекс. size() всегда возвращает количество доступных элементов, начиная с индекса 0 и далее. Когда вы добавляете положительное смещение, size() уменьшается. Добавление отрицательного смещения увеличит size(). Кумулятивное смещение не должно быть ниже нуля, иначе поведение будет неопределенным.

// https://github.com/KubaO/stackoverflown/tree/master/questions/mex-32490874

#include <cstddef>
#include <cstdint>
#include <cmath>
#include <vector>
#include <array>

template <typename T> struct Traits;

template <typename T, class Tr = Traits<T>> class MxArray {
   const mxArray * m_data;
   T * m_ptr;
public:
   typedef T * iterator;
   typedef const T * const_iterator;
   MxArray(const mxArray * data[], int index) :
      m_data(Tr::check(data, index)),
      m_ptr(Tr::ptr(m_data))
   {}
   MxArray(int n, int m, mxArray * out[], index) :
      m_data(out[index] = Tr::create(n, m)),
      m_ptr(Tr::ptr(out[index]))
   {}
   size_t size() const { return mxGetNumberOfElements(m_data) - offset(); }
   int n_dims() const { return mxGetNumberOfDimensions(m_data); }
   int dim(int i) const { return mxGetDimensions(m_data)[i]; }
   inline ptrdiff_t offset() const { return m_ptr - Tr::ptr(m_data); }
   T operator[](ptrdiff_t i) const { return m_ptr[i]; }
   T & operator[](ptrdiff_t i) { return m_ptr[i]; }
   MxArray & operator+=(ptrdiff_t offset) {
      m_ptr += offset;
      return *this;
   }
   friend MxArray operator+(MxArray lhs, ptrdiff_t offset) {
      return lhs += offset;
   }
   void resetOffset() { m_ptr = Tr::ptr(m_data); }
   iterator begin() { return m_ptr; }
   iterator end() { return m_ptr + size(); }
   const_iterator begin() const { return m_ptr; }
   const_iterator end() const { return m_ptr + size(); }
};

Трейты для типов double и mxLogical (не bool!) реализуют проверку типов и передают функциональность соответствующим API-интерфейсам mex.

const char kMsgId[] = "SO:MexExample";

template <> struct Traits<double> {
   static const mxArray * check(const mxArray * data[], int index) {
      if (! mxIsDouble(data[index]))
         mexErrMsgIdAndTxt(kMsgId, "Expected real type for input parameter #%d", index);
      return data[index];
   }
   static double * ptr(const mxArray * data) {
      return mxGetPr(data);
   }
   static mxArray * create(int n, int m) {
      return mxCreateDoubleMatrix(n, m, mxREAL);
   }
};

template <> struct Traits<mxLogical> {
   static const mxArray * check(const mxArray * data[], int index) {
      if (! mxIsLogical(data[index]))
         mexErrMsgIdAndTxt(kMsgId, "Expected logical type for input parameter #%d", index);
      return data[index];
   }
   static mxLogical * ptr(const mxArray * data) {
      return mxGetLogicals(data);
   }
   static mxArray * create(int n, int m) {
      return mxCreateLogicalMatrix(n, m);
   }
};

Наконец, DoubleArray и LogicalArray — это псевдонимы простых типов, которые параметризируют MxArray для нужного типа.

typedef MxArray<double> DoubleArray;
typedef MxArray<mxLogical> LogicalArray;

Неопределенное поведение — обращаться к данным mxLogical, как если бы это были данные bool, но допустимо преобразование между одним значением mxLogical и bool.

Причина возникновения проблем заключалась в том, что не использовался вызов malloc после объявления указателя:

double* Variable;
*Variable = 0.0;
Function_that_Calculates_Variable( Variable );

Это вызывало либо сбои, либо переполнение памяти. Я не понял/забыл, что вам нужно было выделить указатель на одно значение. Поэтому я заменил его на:

double Variable;
Variable = 0.0;
Function_that_Calculates_Variable( &Variable );

который делает то же самое, за исключением того, что не требует вызова malloc и связанного с ним бесплатного, в то же время допуская несколько выходных данных в этой форме из функции.

Задача решена! (После гораздо большего количества часов отладки, чем я готов признать: P)