c - 优化 C 循环

Question

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

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

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

  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”(日晷求解器使用的数据类型)中保存的 vector 的第 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) 最佳做法是尝试将尽可能多的计算移出循环。就像我在上面的循环中所做的一样。