c - 将大方阵乘以转置比大方阵相乘慢...如何解决？

Question

显然，转置一个矩阵然后将其相乘比仅将两个矩阵相乘要快。然而，我的代码现在并没有这样做，我不知道为什么......(正常的乘法只是三重嵌套循环，它给了我大约 1.12 秒来乘以 1000x1000 矩阵，而这段代码给了我 8时间(所以更慢而不是更快)...我现在迷路了任何帮助将不胜感激!:D

A = malloc (size*size * sizeof (double));
B = malloc (size*size * sizeof (double));
C = malloc (size*size * sizeof (double));



/* initialise array elements */
for (row = 0; row < size; row++){
    for (col = 0; col < size; col++){
      A[size * row + col] = rand();
      B[size * row + col] = rand();
    }
  }

t1 = getTime();

/* 这里是要测量的代码 */

T = malloc (size*size * sizeof(double));

for(i = 0; i < size; ++i) {
  for(j = 0; j <= i ; ++j) {
    T[size * i + j] = B[size * j + i];
  }
}

for (j = 0; j < size; ++j) {
  for (k = 0; k < size; ++k) {
    for (m = 0; m < size; ++m) {
      C[size * j + k] = A[size * j + k] * T[size * m + k];
        }
  }
}


t2 = getTime();

Answer 1

我发现了几个问题。

1. 您只是设置C[size * j + k] 的值，而不是递增它。尽管这是计算中的错误，但它不应该影响性能。此外，您需要在最内层循环开始之前将 C[size * j + k] 初始化为 0.0。否则，您将增加一个未初始化的值。这是一个严重的问题，可能会导致溢出。

2. 乘法项错误。

   请记住，您的乘法项需要表示:

         C[j, k] += A[j, m] * B[m, k], which is
         C[j, k] += A[j, m] * T[k, m]
   

   而不是

         C[size * j + k] = A[size * j + k] * T[size * m + k];
   

   你需要

         C[size * j + k] += A[size * j + m] * T[size * k + m];
                     //  ^  ^                 ^^^^^^^^^^^^^^^^
                     //  |  |                 Need to get T[k, m], not T[m, k]
                     //  |  ^^^^^^^^^^^^^^^^
                     //  |  Need to get A[j, m], not A[j, k]
                     //  ^^^^ Increment, not set.
   

我认为，除了错误之外，损害性能的主要罪魁祸首是您对 T[size * m + k] 的使用。当您这样做时，需要大量内存跳转(m 是循环中变化最快的变量)来获取数据。当您使用正确的术语 T[size * k + m] 时，相关内容将会减少，您应该会看到性能改进。

总而言之，使用:

for (j = 0; j < size; ++j) {
   for (k = 0; k < size; ++k) {
      C[size * j + k] = 0.0;
      for (m = 0; m < size; ++m) {
         C[size * j + k] += A[size * j + m] * T[size * k + m];
      }
   }
}

您可以通过使用以下方法获得更高的性能:

double* a = NULL;
double* c = NULL;
double* t = NULL;

for (j = 0; j < size; ++j) {
   a = A + (size*j);
   c = C + (size*j);
   for (k = 0; k < size; ++k) {
      t = T + size*k;
      c[k] = 0.0;
      for (m = 0; m < size; ++m) {
         c[k] += a[m] * t[m];
      }
   }
}

PS我还没有测试过代码。只是给你一些想法。