c++ - 使用 operator() 扩展 Eigen::EigenBase

Question

我想对偏移量为 (0.5, 0) 和 (0, 0.5) 的网格使用 MatrixXd 类。在数学公式中，计算单元格 i,i+1 之间的速度，并将其写为 vel(i+0.5,j)。我想介绍这样的语法:

#include <Eigen/Dense>

int main() {
  Eigen::MatrixXd m = Eigen::MatrixXd::Zero(5,5);
  // Want to use similar syntax:
  // m(0, 1.5) = 1.0;
  // and
  // m(3.5, 1) = 2.0;
  // Instead of:
  m(0, 2) = 1.0;
  m(4, 1) = 2.0;
}

使用 EIGEN_MATRIXBASE_PLUGIN像这个:

inline Scalar& operator()(int r, int c) {
  return Base::operator()(r, c);
}

inline Scalar& operator()(double r, int c) {
  return Base::operator()(int(r + 0.5), c);
}

inline Scalar& operator()(int r, double c) {
  return Base::operator()(r, int(c + 0.5));
}

但是，这种方法:

仅适用于 X 轴或仅 Y 轴偏移，不能同时用于两者。

仅适用于硬编码到插件中的特定偏移量。

打破一些内部特征对流，可以是demostrated通过尝试使用 IncompleteLUT 预处理器编译 BiCG 示例:

  int n = 10000;
  VectorXd x(n), b(n);
  SparseMatrix<double> A(n,n);
  /* ... fill A and b ... */ 
  BiCGSTAB<SparseMatrix<double>,IncompleteLUT<double>> solver;
  solver.compute(A);
  x = solver.solve(b);

导致以下错误:

term does not evaluate to a function taking 1 arguments
'Eigen::SparseMatrix<double,1,int>::insertBackByOuterInnerUnordered': function does not take 1 arguments

添加 operator()(double offset_col, double offset_row) 来解决第二个问题，如下所示:

double r_offset = -0.5, c_offset = -0.5;

inline void set_r_offset(double val) { r_offset = val; }
inline void set_c_offset(double val) { c_offset = val; }

inline double get_r_offset() { return r_offset; }
inline double get_c_offset() { return c_offset; }

inline Scalar& operator()(double r, double c) {
//  double r_offset = -0.5, c_offset = -0.5;
  return Base::operator()(int(r - r_offset), int(c - c_offset));
}

这导致非法免费:

==6035== Invalid free() / delete / delete[] / realloc()
==6035==    at 0x4C30D3B: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==6035==    by 0x4E4224A: aligned_free (Memory.h:177)
==6035==    by 0x4E4224A: conditional_aligned_free<true> (Memory.h:230)
==6035==    by 0x4E4224A: conditional_aligned_delete_auto<double, true> (Memory.h:416)
==6035==    by 0x4E4224A: resize (DenseStorage.h:406)
==6035==    by 0x4E4224A: resize (PlainObjectBase.h:293)
==6035==    by 0x4E4224A: resize_if_allowed<Eigen::Matrix<double, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> >, double, double> (AssignEvaluator.h:720)
==6035==    by 0x4E4224A: call_dense_assignment_loop<Eigen::Matrix<double, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> >, Eigen::internal::assign_op<double, double> > (AssignEvaluator.h:734)
==6035==    by 0x4E4224A: run (AssignEvaluator.h:879)
==6035==    by 0x4E4224A: call_assignment_no_alias<Eigen::Matrix<double, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> >, Eigen::internal::assign_op<double, double> > (AssignEvaluator.h:836)
==6035==    by 0x4E4224A: call_assignment<Eigen::Matrix<double, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> >, Eigen::internal::assign_op<double, double> > (AssignEvaluator.h:804)
==6035==    by 0x4E4224A: call_assignment<Eigen::Matrix<double, -1, -1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> > > (AssignEvaluator.h:782)
==6035==    by 0x4E4224A: _set<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> > > (PlainObjectBase.h:710)
==6035==    by 0x4E4224A: operator=<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, -1> > > (Matrix.h:225)
==6035==    by 0x11044C: main (Runner.cpp:16)
==6035==  Address 0x2e642f73726573 is not stack'd, malloc'd or (recently) free'd

如果偏移量不是作为类成员引入的，而是 operator() 中的局部变量，则 valgrind 不会检测到错误。

是否可以使用可设置的偏移量实现新的 MatrixXd::operator()(double, double)？

编辑:
Operator() 在父类 DenseCoeffsBase 中定义:

EIGEN_DEVICE_FUNC
EIGEN_STRONG_INLINE CoeffReturnType operator()(Index row, Index col) const
{
   eigen_assert(row >= 0 && row < rows()
       && col >= 0 && col < cols());
   return coeff(row, col);
}

Answer 1

也许，我发现您的运算符存在一个问题，它返回对 Scalar 临时对象的引用:

inline Scalar& operator()(double r, double c) {
    //  double r_offset = -0.5, c_offset = -0.5;
    return Base::operator()(int(r - r_offset), int(c - c_offset));
}

因此，您应该通过拷贝返回 Scalar。

你能分享 Base::operator()(int(r - r_offset), int(c - c_offset)); 的代码吗？