Блог

Я пытаюсь реализовать код, показанный в этом pdf. Более точно (стр. 50):

 #define SM (CLS / sizeof (double))

for (i = 0; i < N; i  = SM)
  for (j = 0; j < N; j  = SM)
    for (k = 0; k < N; k  = SM)
      for (i2 = 0, rres = amp;res[i][j],
           rmul1 = amp;mul1[i][k]; i2 < SM;
             i2, rres  = N, rmul1  = N)
         for (k2 = 0, rmul2 = amp;mul2[k][j];
              k2 < SM;   k2, rmul2  = N)
           for (j2 = 0; j2 < SM;   j2)
              rres[j2]  = rmul1[k2] * rmul2[j2];
 

Теперь, насколько я понимаю, rres есть int* , rmul2 и rmul1 тоже.

Должно ли это выглядеть так

 int *rres;
int *rmul1;
int *rmul2;

for (i = 0; i < N; i  = SM)
   for (j = 0; j < N; j  = SM)
      for (k = 0; k < N; k  = SM)
         for (i2 = 0, *rres = amp;res[i][j],
              *rmul1 = amp;mul1[i][k]; i2 < SM;
                i2, rres  = N, rmul1  = N)
            for (k2 = 0, *rmul2 = amp;mul2[k][j];
                 k2 < SM;   k2, rmul2  = N)
               for (j2 = 0; j2 < SM;   j2)
                  rres[j2]  = rmul1[k2] * rmul2[j2];
 

Поскольку это кажется мне более или менее разумным, бит дает неправильные результаты. Например, если у меня есть две матрицы 2 x 2, случайные значения 0 или 1, я получаю:

 -1520010527 23996350 
212687419 207125308
 

что далеко не хорошо. Я предполагаю, что я использую * неправильно, но я не могу сказать, где…

РЕДАКТИРОВАТЬ: объявление res :

   int **res = (int **)malloc(N * sizeof(int *));
  for (int i = 0; i < N; i  ) {
      res[i] = (int *)malloc(N * sizeof(int));
  }
  for (int i = 0; i < N; i  ) {
      for (int j = 0; j < N; j  ) {
          res[i][j] = 0;
      }
  }
 

Ваше предположение о типе rres rmul2 и rmul1 кажется неправильным: элементы матрицы, вероятно, имеют double тип, а выражения в начальных предложениях for циклов определенно не должны иметь инициала * .

Тем не менее, этот код умножения матрицы должен работать одинаково для любого типа элемента, но оптимизация кэша, объясненная Ульрихом Дреппером, делает другое предположение: N должно быть кратным SM , и это вряд ли будет иметь место с #define SM (CLS / sizeof(double)) и N = 2 в вашем тестовом примере, что может объяснить ваши наблюдения.

Код должен быть:

 #define SM  (CLS / sizeof(double))

void matrix_multiply(double res[N][N], const double mul1[N][N], const double mul2[N][N]) {
    size_t i, j, k, i2, j2, k2;
    double *rres, *rmul1, *rmul2;
    
    assert(N % SM == 0);

    for (i = 0; i < N; i  = SM) {
        for (j = 0; j < N; j  = SM) {
            for (k = 0; k < N; k  = SM) {
                for (i2 = 0, rres = amp;res[i][j], rmul1 = amp;mul1[i][k];
                     i2 < SM;   i2, rres  = N, rmul1  = N) {
                    for (k2 = 0, rmul2 = amp;mul2[k][j];
                         k2 < SM;   k2, rmul2  = N) {
                        for (j2 = 0; j2 < SM;   j2)
                            rres[j2]  = rmul1[k2] * rmul2[j2];
                    }
                }
            }
        }
    }
}