0

令 A 为对称矩阵，令 v 为向量。我从 A 中提取从 j 开始的 n 列块并将其乘以 v 使用

VectorXd a;
a = A.middleCols(j,n).selfadjointView<Lower>() * v // does not compile

因为这不能编译，而这

a = MatrixXd(A.middleCols(j,n).selfadjointView<Lower>()) * v

确实，我想知道第二个版本是否复制了

A.middleCols(j,n).selfadjointView<Lower>()

还是直接执行计算？

感谢您的任何提示。

编辑：我怀疑这个问题与参数类型有关，因为我得到了错误：

invalid argument type 'typename ConstSelfAdjointViewReturnType.... to unary expression'

实际上， A 是使用 const 引用传递的函数的参数

const MatrixXd& A 
const Ref<const MatrixXd>& A

这是一个例子：

// this version doesn't compile
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = A.middleCols(3, 4).selfadjointView<Lower>() * v;
MatrixXd M(code_fill(){...}); 
// code_fill is the function filling the lower triangular part of a symmetric matrix
M.block(1, 2, 3, 4).triangularView<Lower>() += B.selfadjointView<Lower>();

return M;
}

// this version compiles
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = MatrixXd(A.middleCols(3, 4).selfadjointView<Lower>()) * v;
MatrixXd M(code_fill(){...}); 
// code_fill is the function filling the lower triangular part of a symmetric matrix
Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

return M;
}

EDIT2关于我最初的问题和我在编辑部分添加的示例，我对复制有点困惑。据我了解工作版本和非工作版本之间的区别，该行

Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

之所以有效，是因为它的 lhs 告诉 Eigen M.block(1, 2, 3, 4).triangularView() 实际上是一个矩阵，而不是对矩阵的引用。否则，运算符 += 将通过一个错误，即该运算符没有为 .block() 重载。所以我最初的问题是 Matrix(...) 是否只告诉它是一个 Matrix 来启用计算，或者更确切地说将 ... 复制到一个 Matrix 中？谢谢！

4

1 回答 1

1

以下表达式：

A.middleCols(j,n).selfadjointView<Lower>()

不创建任何副本。

另一方面，为了避免产品结果的临时性，您可以添加.noalias()：

a.noalias() = M.middleCols(j,n).selfadjointView<Lower>() * v;

这只需要立即分配密集产品以允许以下代码：

a = M * a;

按预期行事。

编辑：

关于您的编译问题，以下编译良好：

#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  int n = 10;
  MatrixXd M = MatrixXd::Random(n,2*n);
  VectorXd v = VectorXd::Random(n);
  VectorXd a;
  a.noalias() = M.middleCols(2,n).selfadjointView<Upper>() * v;

  const Ref<const MatrixXd>& A = M;
  a.noalias() = A.middleCols(2,n).selfadjointView<Upper>() * v;
}

编辑2

以下行：

MatrixXd(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

没有意义，因为您正在创建要分配的临时副本。回想一下，这里MatrixXd(whatever)调用了Eigen::Matrix构造函数。将自伴随视图分配给三角形视图也没有任何意义。我什至想不出什么是合理的行为。如果您只想更新下三角形部分，请执行以下操作：

M.block(1, 2, 3, 4).triangularView<Lower>() += B;

在这种情况下，triangularView充当书写掩码。

于 2017-01-30T13:09:01.047 回答