Расширение Eigen :: EigenBase с помощью оператора () - PullRequest
Новая
       11

Расширение Eigen :: EigenBase с помощью оператора ()

0 голосов
/ 02 февраля 2020

Я хотел бы использовать класс 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. Работает только для указанного смещения c, жестко закодированного в плагин.
  3. Breaks некоторые внутренние собственные конвекции, которые можно демострировать , пытаясь скомпилировать пример BiCG с прекондиционером IncompleteLUT: 
  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

Добавление оператора () (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) с настраиваемыми смещениями?

EDIT : Оператор () определен в родительском классе 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);
}

1 Ответ

0 голосов
/ 02 февраля 2020

Возможно, я вижу одну проблему с вашим оператором, который возвращает ссылку на временный объект 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));?

...