Оптимизация MEX-файла для ускорения умножения больших многомерных массивов (узкое место кода)

Question

Я пытаюсь выполнить такие операции, как умножение массивов 7D, которые содержат 32 миллиона элементов. Я написал файл MEX, так как у меня сложилось впечатление, что эти операции должны выполняться быстрее в C, чем в Matlab. Тем не менее, я обнаружил, что MEX-файл примерно в два раза медленнее, чем выполнение операций непосредственно в Matlab (2017b).

Пример операции, которую я хочу выполнить:

T8  = rand(1,1e3,2,2,2,2,2);
wsm = rand(1e3,1e3,2,2);
CM  = bsxfun(@times,T8,wsm);

На моей машине это занимает 0.117065 секунд (я называю это и другие подобные операции ~1000 раз за прогон модели, а модель запускается тысячи раз для оптимизации параметров - эти операции делают оптимизацию слишком медленной).

Вот MEX-файл, который я написал, он использует 7 для циклов для доступа к элементам T8 и wsm посредством линейного индексирования (может быть, я должен был бы получить доступ к элементам более эффективным способом или избегать циклов?):

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize i, j, k, l, m, n, o, I, J, K, L, M, N, O;
    mwSize *dims,*dims1;
    double *T8, *wsm, *CM;
      T8  = mxGetPr(prhs[0]);
      wsm = mxGetPr(prhs[1]);

      dims = mxGetDimensions(prhs[0]);
      dims1 = mxGetDimensions(prhs[1]);
      dims[0] = dims1[0];

      I = dims[0];
      J = dims[1];
      K = dims[2];
      L = dims[3];
      M = dims[4];
      N = dims[5];
      O = dims[6];

      plhs[0] = mxCreateNumericArray(7,dims,mxDOUBLE_CLASS,mxREAL);
      CM = mxGetPr(plhs[0]);

      for( o=0; o<O; o++ ) {
          for( n=0; n<N; n++ ) {
              for( m=0; m<M; m++ ) {
                  for( l=0; l<L; l++ ) {
                      for( k=0; k<K; k++ ) {
                          for( j=0; j<J; j++ ) {
                              for( i=0; i<I; i++ ) {
                                  *CM++ = T8[j + k*J + +l*J*K + m*L*J*K + n*M*L*J*K + o*N*M*L*J*K] * wsm[i + j*I + k*I*J + l*I*J*K];
                              }
                          }
                      }
                  }
              }
          }
      }
}

Когда я вызываю вышеупомянутый файл MEX

CM = arrayProduct(T8,wsm);

это занимает 0,215211 секунды (почти вдвое дольше).

Мой код был очень свободно основан на коде, предложенном здесь ( https://uk.mathworks.com/matlabcentral/answers/210352-optimize-speed-up-a-big-and-slow-matrix-operation-with-addition-and-bsxfun).

Будем очень благодарны за любые предложения относительно того, что я могу сделать по-другому для ускорения работы моего кода!

Answer 1

Это большая ошибка, предполагать, что вы можете победить Matlab в такой простой математике. Matlab оптимизирован с самого начала для выполнения математических матриц.

Иногда есть веские причины для написания MEX-функций, в том числе из соображений производительности, но обычно это происходит в тех случаях, когда чистое решение Matlab невозможно написать оптимальным способом (например, когда вам нужно написать множество явных циклов).

Две основные причины, по которым ваш код может быть медленнее, чем математика оптимизированной матрицы, уже присутствующая в Matlab:

1. Matlab может использовать несколько потоков для параллельных вычислений. Ваш код не, но действительно оптимальное решение, вероятно, будет.
2. Возможно, вы допустили ошибку в схеме доступа к памяти, что привело к снижению частоты обращений к кешу.

Еще один способ взглянуть на это так: если Matlab нельзя доверять для реализации умножения оптимальным образом, будут ли люди использовать его для серьезной математики на больших наборах данных? Есть алгоритмы, которых Matlab не знает, и иногда их можно ускорить с помощью MEX, но умножение не является одним из них.