使用 operator(） 扩展 Eigen::EigenBase

我想将MatrixXd类用于偏移量为(0.5,0(和(0,0.5(的网格。在数学公式中，速度是在单元格 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));
}

但是，这种方法：

1. 仅适用于 X 轴或 Y 轴偏移，不能同时用于两者。
2. 仅适用于硬编码到插件中的特定偏移量。
3. 打破了一些内部特征对流，可以通过尝试使用 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

添加运算符(((双offset_col，双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);
}