Матричное умножение в файле cpp для Matlab

Question

Как бы сделать матричное умножение в формате cpp, который будет скомпилирован в файл mex?

Мое нормальная матрица умножения в сценарии Matlab выглядит следующим образом:

cMatrix = (1/r) * pfMatrix * wcMatrix; %here pfMatrix is 2x3 and wcMatrix is 3x8 
% Hence cMatrix is 2x8 
% r is a scalar 

pfMatrix, wcMatrix и г правильно объявлено в файле cpp, и они имеют то же значение, что и в сценарии. Однако cMatrix не дает мне таких же результатов. Здесь реализация матричного умножения в cpp:

int i, n, j; 
    for (i = 0; i<1; i++) 
    { 
     for (n = 0; n<7; n++) 
     { 
      for (j = 0; j<2; j++) 
      { 
       d->cMatrix[i][n] += (d->pfMatrix[i][j]) * (d->wcMatrix[j][n]); 
      } 
      d->cMatrix[i][n] = (1/d->r) * d->cMatrix[i][n]; 
     } 
    } 

Edit:

Я изменил цикл следующий Ben Voigt ответ. Результаты в cMatrix по-прежнему не идентичны результатам, рассчитанным по сценарию Matlab.

Например:

pfMatrix = [7937.91049469652,0,512;0,7933.81033431703,384]; 
wcMatrix = [-0.880633810389421,-1.04063381038942,-1.04063381038942,-0.880633810389421,-0.815633810389421,-1.10563381038942,-1.10563381038942,-0.815633810389421;-0.125,-0.125,0.125,0.125,-0.29,-0.29,0.29,0.29;100,100,100,100,100,100,100,100]; 
r = 100; 

В этом случае cMatrix(1,1) является:

(pfMatrix(1,1)*wcMatrix(1,1) + pfMatrix(1,2)*wcMatrix(2,1) + pfMatrix(1,3)*wcMatrix(3,1))/r = 442.09 

Однако с mex файла эквивалентный результат является 959.

Edit # 2:

Я нашел ошибку в элементе pfMatrix, который не был правильно объявлены (отсутствует деление на 2). Поэтому ответ Ben Voigt работает правильно. Однако между этими двумя результатами все еще есть небольшая разница (скрипт Matlab дает 442, а mex дает 447, может ли это быть результатом разного типа данных?).

Edit # 3:

Нашли ошибку, и это не было связано с циклом умножения матриц.

Answer 1

Использование вашей матрицы результатов в качестве пространства для царапин - отличная идея. Компилятор должен беспокоиться об aliasing, что означает, что он не может оптимизировать.

Попробуйте явную рабочую переменную, которая также обеспечивает удобное место для обнулить его:

for (int i = 0; i < 2; ++i) { 
    for (int n = 0; n < 8; ++n) { 
     double accum = 0.0; 
     for (int j = 0; j < 3; ++j) { 
      accum += (d->pfMatrix[i][j]) * (d->wcMatrix[j][n]); 
     } 
     d->cMatrix[i][n] = accum/d->r; 
    } 
} 

Ваши диапазоны были также неправильно, что я фиксированы.

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

Answer 2

Умножение между матрицами должно быть таким образом, : A [m] [n] * B [n] [p] = R [m] [p].

Условия, которые вы написали в цикле for, неверны и не учитывают размеры матрицы.

Посмотрите также на собственные библиотеки, которые являются с открытым исходным кодом и обеспечивают простой способ выполнения матричных умножений.