c++ - 矩阵 vector 乘法与 Eigen 系数乘积的结合

Question

我正在尝试使用 Eigen 中的函数来压缩一些代码，特别是那些允许您对数组进行系数乘积的代码，但我可能没有正确使用它们。基本思想包含在这个扩展的 for 循环中:

  for (int dir = 0;  dir < NDIM; dir++)
  {
    for (int i = 0; i < nlocal; i++)
    {
      for (int qp = 0; qp < nVolQuad; qp++)
        qNewPtr[i] += 
          volQuad.weights(qp)*bigStoredVolMatrices[dir](i,qp)*alpha(dir,qp)*fAtQuad(qp);
    }
  }

我想浓缩成:

  for (int dir = 0;  dir < NDIM; dir++)
  {  
    resultVectorDir[dir].noalias() = bigStoredVolMatrices[dir]*
      (volQuad.weights.array()*fAtQuad.array()*alpha.row(dir).array()).matrix();  
  }

  for (int i = 0; i < nlocal; i++)
  {
    for (int dir = 0;  dir < NDIM; dir++)
      qNewPtr[i] += resultVectorDir[dir](i);
  }

或者不在数组和矩阵之间切换，使用类似的东西:

  for (int dir = 0;  dir < NDIM; dir++)
  {  
    resultVectorDir[dir].noalias() = bigStoredVolMatrices[dir]*
      (volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir))));  
  }

  for (int i = 0; i < nlocal; i++)
  {
    for (int dir = 0;  dir < NDIM; dir++)
      qNewPtr[i] += resultVectorDir[dir](i);
  }

对我来说奇怪的是有时这行得通。代码有时会产生所需的输出，但有时也会产生 NaN。我想我可能需要明确地将压缩版本中的 resultVectorDir 归零，但这并没有解决问题。我认为执行此操作顺序可能只是一些微妙的事情？如果能提供任何帮助，我们将不胜感激。

作为这个问题的附录，我用老式的打印语句解决了这个问题，发现我一定没有正确使用数组函数的系数乘积。例如，在 nVolQuad = 9 的情况下，我运行了这段代码:

  for (int dir = 0;  dir < NDIM; dir++)
  {
    tempArray = volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir)));
    for (int qp = 0; qp < nVolQuad; qp++)
    {
      std::cout << std::endl;
      std::cout << tempArray(qp) << " "; 
      std::cout << volQuad.weights(qp)*alpha(dir,qp)*fAtQuad(qp) << " ";
      std::cout << std::endl;
    }
    std::cout << std::endl;
    std::cout << std::endl;
  }

一部分输出如下所示:

-2.23064e-05 -2.23064e-05

1.49458e-154 -3.56902e-05

6.94729e-310 -2.23064e-05

-2.68156e+154 -2.0672e-05

6.94729e-310 -3.30752e-05

6.94729e-310 -2.0672e-05

2.13645e-314 -2.99114e-06

6.94729e-310 -4.78582e-06

6.94729e-310 -2.99114e-06

输出的其他部分类似。第一个条目是正确的，但 tempArray 的其他 8 个条目是无意义的。 tempArray 在循环之前被初始化为 0.0，所以我有点不知所措。我将继续深入研究 Eigen 的文档，以确保我在使用此功能时没有做一些极其愚蠢的事情。

Answer 1

关键错误在于假设当扩展循环正在访问不同大小的数据时，您可以立即从扩展循环转到系数乘积。这里:

volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir)))

由于 volQuad.weights 和 fAtQuad 是使用 Eigen::VectorXd 初始化的，因此它们是列 vector ，但是通过使用 alpha.row(dir)，该特定数据结构是行 vector 。因此，系数明智的产品没有意义，你只会让第一个条目正确。可以通过将语法更改为以下方式轻松解决此问题:

volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir).transpose()))