使用operator（）扩展Eigen :: EigenBase

Question

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

int j;
for (double i=0.5; i<nx; i+=1.0){
   //calculate velocity
   double val = f(h(i-0.5,j), h(i+0.5,j));
   vel(i,j) = val;
}

像这样添加运算符（）（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));
}

所以您应该通过副本退还标量。

您能否共享Base :: operator（）（int（r-r_offset），int（c-c_offset））;？