优化 C 循环

Question

我是多年使用 Matlab 进行数值编程的 C 新手。我已经开发了一个程序来解决一个大型的微分方程系统，但我很确定我做了一些愚蠢的事情，因为在分析代码之后，我惊讶地发现三个循环占用了大约 90% 的计算时间，尽管他们正在执行程序中最琐碎的步骤。

基于这些昂贵的循环，我的问题分为三个部分：

• 将数组初始化为零。当 J 被声明为双精度数组时，数组的值是否初始化为零？如果没有，有没有一种快速的方法将所有元素设置为零？

  void spam(){
      double J[151][151];    
      /* Other relevant variables declared */
      calcJac(data,J,y);
      /* Use J */
  }
  
  static void calcJac(UserData data, double J[151][151],N_Vector y)
  {
      /* The first expensive loop */
      int iter, jter;
      for (iter=0; iter<151; iter++) {
          for (jter = 0; jter<151; jter++) {
              J[iter][jter] = 0;
          }
      }
     /* More code to populate J from data and y that runs very quickly */
  }
  

• 在求解过程中，我需要求解由 P = I - gamma*J 定义的矩阵方程。 P 的构造比求解它定义的方程组花费的时间更长，所以我正在做的事情可能是错误的。在下面相对较慢的循环中，访问包含在结构“数据”中的矩阵是慢速组件还是与循环有关的其他内容？

  for (iter = 1; iter<151; iter++) {
      for(jter = 1; jter<151; jter++){
          P[iter-1][jter-1] = - gamma*(data->J[iter][jter]);
      }
  }
  

• 是否有矩阵乘法的最佳实践？在下面的循环中，Ith(v,iter) 是一个宏，用于获取保存在 N_Vector 结构“v”（日晷求解器使用的数据类型）中的向量的第 iter 个分量。特别是，有没有最好的方法来获得 v 和 J 的行之间的点积？

  Jv_scratch = 0;
  int iter, jter;
  for (iter=1; iter<151; iter++) {
      for (jter=1; jter<151; jter++) {
          Jv_scratch += J[iter][jter]*Ith(v,jter);
      }
      Ith(Jv,iter) = Jv_scratch;
      Jv_scratch = 0;
  }
  

Answer 1

1) 不，他们不是你可以按如下方式 memset 数组：

memset( J, 0, sizeof( double ) * 151 * 151 );

或者你可以使用数组初始化器：

double J[151][151] = { 0.0 };

2) 好吧，您正在使用相当复杂的计算来计算 P 的位置和 J 的位置。

您可能会获得更好的性能。通过作为指针单步执行：

for (iter = 1; iter<151; iter++) 
{
    double* pP = (P - 1) + (151 * iter);
    double* pJ = data->J + (151 * iter);

    for(jter = 1; jter<151; jter++, pP++, pJ++ )
    {
         *pP = - gamma * *pJ;
    }
}

通过这种方式，您可以将各种数组索引计算移到循环之外。

3) 最佳做法是尝试将尽可能多的计算移出循环。就像我在上面的循环中所做的一样。

Answer 2

首先，我建议您将问题分成三个单独的问题。很难回答所有三个问题。例如，我对数值分析的工作不多，所以我只回答第一个。

首先，堆栈上的变量没有为您初始化。但是有更快的方法来初始化它们。在你的情况下，我建议使用 memset：

static void calcJac(UserData data, double J[151][151],N_Vector y)
{
   memset((void*)J, 0, sizeof(double) * 151 * 151);
   /* More code to populate J from data and y that runs very quickly */
}

memset 是一个快速的库例程，用于用特定模式的字节填充内存区域。碰巧将double 的所有字节设置为零会将double 设置为零，因此请利用库的快速例程（可能会用汇编程序编写以利用SSE 之类的东西）。

Answer 3

其他人已经回答了您的一些问题。关于矩阵乘法的主题；很难为此编写一个快速的算法，除非您对缓存体系结构等非常了解（速度慢是由于您访问数组元素的顺序会导致数千次缓存未命中）。

如果您想了解“matrix-multiplication”、“cache”、“blocking”等术语，可以尝试在 Google 上搜索快速库中使用的技术。但我的建议是，如果性能是关键，就使用预先存在的数学库。

Answer 4

将数组初始化为零。 当 J 被声明为双精度 数组是数组的值 初始化为零？如果没有，有没有 将所有元素设置为的快速方法 零？

这取决于数组的分配位置。如果它在文件范围内声明或声明为静态，则 C 标准保证所有元素都设置为零。如果您在初始化时将第一个元素设置为一个值，则可以保证相同，即：

double J[151][151] = {0}; /* set first element to zero */

通过将第一个元素设置为某个值，C 标准保证数组中的所有其他元素都设置为零，就好像数组是静态分配的一样。

实际上对于这种特定情况，我非常怀疑无论您使用哪个系统，在堆栈上分配 151*151*sizeof(double) 字节是否明智。您可能必须动态分配它，然后以上都不重要。然后您必须使用 memset() 将所有字节设置为零。

在 下面比较慢的循环，是 访问包含的矩阵 在结构“数据”中缓慢 组件还是其他东西 关于循环？

您应该确保从中调用的函数是内联的。否则，您无法优化循环：什么是最佳的高度依赖于系统（即如何构建物理高速缓存）。最好将这种优化留给编译器。

您当然可以通过手动优化来混淆代码，例如倒数而不是倒数，或者使用 ++i 而不是 i++ 等。但是编译器确实应该能够为您处理这些事情。

至于矩阵加法，我不知道数学上最有效的方法，但我怀疑它与代码效率的关系不大。这里的大盗是双重类型。除非你真的需要高精度，否则我会考虑使用 float 或 int 来加速算法。