矩阵-向量乘法与特征系数积的结合

Combination of matrix-vector multiplication with coefficient-wise product in eigen

本文关键字：结合 特征 向量 矩阵      更新时间：2023-10-16

我正试图使用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);
  }

对我来说奇怪的是有时这是有效的。代码有时会产生期望的输出，但有时也会产生纳米网络。我认为我可能需要显式地将浓缩版本中的resultVectorDir归零，但这并没有解决问题。我想执行这个操作顺序可能有些微妙的地方?如能提供任何帮助，我将不胜感激。

作为这个问题的补充，我用老式的print语句解决了这个问题，发现我一定没有正确地使用数组函数的系数明智乘积。例如，在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.23064 e-05 -2.23064 e-05

1.49458 e - 154 -3.56902 e-05

6.94729 e - 310 -2.23064 e-05

-2.68156 e + 154 -2.0672 e-05

6.94729 e - 310 -3.30752 e-05

6.94729 e - 310 -2.0672 e-05

2.13645 e - 314 -2.99114 e-06

6.94729 e - 310 -4.78582 e-06

6.94729 e - 310 -2.99114 e-06

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

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

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

。weights和fAtQuad是用Eigen::VectorXd初始化的，它们是列向量，但是通过使用alpha.row(dir)，这个特定的数据结构是行向量。因此，系数和乘积是没有意义的你只会得到第一个正确的项。通过将语法更改为:

，可以很容易地解决这个问题。
volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir).transpose()))
相关文章：