From 8657f3083a633d6d4a74e9bcecd3856dd3b513b6 Mon Sep 17 00:00:00 2001 From: Jason Rhinelander Date: Thu, 4 Aug 2016 13:21:39 -0400 Subject: [PATCH] Fix eigen copying of non-standard stride values Some Eigen objects, such as those returned by matrix.diagonal() and matrix.block() have non-standard stride values because they are basically just maps onto the underlying matrix without copying it (for example, the primary diagonal of a 3x3 matrix is a vector-like object with .src equal to the full matrix data, but with stride 4). Returning such an object from a pybind11 method breaks, however, because pybind11 assumes vectors have stride 1, and that matrices have strides equal to the number of rows/columns or 1 (depending on whether the matrix is stored column-major or row-major). This commit fixes the issue by making pybind11 use Eigen's stride methods when copying the data. --- example/eigen.cpp | 10 ++++++++++ example/eigen.py | 10 ++++++++++ example/eigen.ref | 17 +++++++++++++++++ include/pybind11/eigen.h | 10 +++++----- 4 files changed, 42 insertions(+), 5 deletions(-) diff --git a/example/eigen.cpp b/example/eigen.cpp index 728b575fd..9c33c551e 100644 --- a/example/eigen.cpp +++ b/example/eigen.cpp @@ -56,6 +56,16 @@ void init_eigen(py::module &m) { m.def("cholesky5", &cholesky5); m.def("cholesky6", &cholesky6); + // Returns diagonals: a vector-like object with an inner stride != 1 + m.def("diagonal", [](const Eigen::Ref &x) { return x.diagonal(); }); + m.def("diagonal_1", [](const Eigen::Ref &x) { return x.diagonal<1>(); }); + m.def("diagonal_n", [](const Eigen::Ref &x, int index) { return x.diagonal(index); }); + + // Return a block of a matrix (gives non-standard strides) + m.def("block", [](const Eigen::Ref &x, int start_row, int start_col, int block_rows, int block_cols) { + return x.block(start_row, start_col, block_rows, block_cols); + }); + m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); }); diff --git a/example/eigen.py b/example/eigen.py index 6cdc3940b..04078b9d9 100644 --- a/example/eigen.py +++ b/example/eigen.py @@ -12,6 +12,8 @@ from example import sparse_passthrough_r, sparse_passthrough_c from example import double_row, double_col from example import double_mat_cm, double_mat_rm from example import cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6 +from example import diagonal, diagonal_1, diagonal_n +from example import block try: import numpy as np import scipy @@ -78,3 +80,11 @@ for chol in [cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6]: print("cholesky" + str(i) + " " + ("OK" if (mymat == np.array([[1,0,0], [2,3,0], [4,5,6]])).all() else "NOT OKAY")) i += 1 +print("diagonal() %s" % ("OK" if (diagonal(ref) == ref.diagonal()).all() else "FAILED")) +print("diagonal_1() %s" % ("OK" if (diagonal_1(ref) == ref.diagonal(1)).all() else "FAILED")) +for i in range(-5, 7): + print("diagonal_n(%d) %s" % (i, "OK" if (diagonal_n(ref, i) == ref.diagonal(i)).all() else "FAILED")) + +print("block(2,1,3,3) %s" % ("OK" if (block(ref, 2, 1, 3, 3) == ref[2:5, 1:4]).all() else "FAILED")) +print("block(1,4,4,2) %s" % ("OK" if (block(ref, 1, 4, 4, 2) == ref[1:, 4:]).all() else "FAILED")) +print("block(1,4,3,2) %s" % ("OK" if (block(ref, 1, 4, 3, 2) == ref[1:4, 4:]).all() else "FAILED")) diff --git a/example/eigen.ref b/example/eigen.ref index 93e88adb9..8ccd1f472 100644 --- a/example/eigen.ref +++ b/example/eigen.ref @@ -33,3 +33,20 @@ cholesky3 OK cholesky4 OK cholesky5 OK cholesky6 OK +diagonal() OK +diagonal_1() OK +diagonal_n(-5) OK +diagonal_n(-4) OK +diagonal_n(-3) OK +diagonal_n(-2) OK +diagonal_n(-1) OK +diagonal_n(0) OK +diagonal_n(1) OK +diagonal_n(2) OK +diagonal_n(3) OK +diagonal_n(4) OK +diagonal_n(5) OK +diagonal_n(6) OK +block(2,1,3,3) OK +block(1,4,4,2) OK +block(1,4,3,2) OK diff --git a/include/pybind11/eigen.h b/include/pybind11/eigen.h index b9f22b049..16ccaa845 100644 --- a/include/pybind11/eigen.h +++ b/include/pybind11/eigen.h @@ -40,8 +40,8 @@ public: static constexpr bool value = decltype(test(std::declval()))::value; }; -// Eigen::Ref satisfies is_eigen_dense, but isn't constructible, which means we can't load -// it (since there is no reference!), but we can cast from it. +// Eigen::Ref satisfies is_eigen_dense, but isn't constructible, so it needs a special +// type_caster to handle argument copying/forwarding. template class is_eigen_ref { private: template static typename std::enable_if< @@ -126,7 +126,7 @@ struct type_caster::value && /* Buffer dimensions */ { (size_t) src.size() }, /* Strides (in bytes) for each index */ - { sizeof(Scalar) } + { sizeof(Scalar) * src.innerStride() } )).release(); } else { return array(buffer_info( @@ -142,8 +142,8 @@ struct type_caster::value && { (size_t) src.rows(), (size_t) src.cols() }, /* Strides (in bytes) for each index */ - { sizeof(Scalar) * (rowMajor ? (size_t) src.cols() : 1), - sizeof(Scalar) * (rowMajor ? 1 : (size_t) src.rows()) } + { sizeof(Scalar) * src.rowStride(), + sizeof(Scalar) * src.colStride() } )).release(); } }