Eigen<->numpy referencing support

This commit largely rewrites the Eigen dense matrix support to avoid
copying in many cases: Eigen arguments can now reference numpy data, and
numpy objects can now reference Eigen data (given compatible types).

Eigen::Ref<...> arguments now also make use of the new `convert`
argument use (added in PR #634) to avoid conversion, allowing
`py::arg().noconvert()` to be used when binding a function to prohibit
copying when invoking the function.  Respecting `convert` also means
Eigen overloads that avoid copying will be preferred during overload
resolution to ones that require copying.

This commit also rewrites the Eigen documentation and test suite to
explain and test the new capabilities.

docs/advanced/cast/eigen.rst

@@ -1,48 +1,308 @@
 Eigen
-=====
+#####
 
 `Eigen <http://eigen.tuxfamily.org>`_ is C++ header-based library for dense and
 sparse linear algebra. Due to its popularity and widespread adoption, pybind11
-provides transparent conversion support between Eigen and Scientific Python linear
+provides transparent conversion and limited mapping support between Eigen and
-algebra data types.
+Scientific Python linear algebra data types.
 
-Specifically, when including the optional header file :file:`pybind11/eigen.h`,
+To enable the built-in Eigen support you must include the optional header file
-pybind11 will automatically and transparently convert
+:file:`pybind11/eigen.h`.
 
-1. Static and dynamic Eigen dense vectors and matrices to instances of
+Pass-by-value
-   ``numpy.ndarray`` (and vice versa).
+=============
 
-2. Returned matrix expressions such as blocks (including columns or rows) and
+When binding a function with ordinary Eigen dense object arguments (for
-   diagonals will be converted to ``numpy.ndarray`` of the expression
+example, ``Eigen::MatrixXd``), pybind11 will accept any input value that is
-   values.
+already (or convertible to) a ``numpy.ndarray`` with dimensions compatible with
+the Eigen type, copy its values into a temporary Eigen variable of the
+appropriate type, then call the function with this temporary variable.
 
-3. Returned matrix-like objects such as Eigen::DiagonalMatrix or
+Sparse matrices are similarly copied to or from
-   Eigen::SelfAdjointView will be converted to ``numpy.ndarray`` containing the
+``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` objects.
-   expressed value.
 
-4. Eigen sparse vectors and matrices to instances of
+Pass-by-reference
-   ``scipy.sparse.csr_matrix``/``scipy.sparse.csc_matrix`` (and vice versa).
+=================
 
-This makes it possible to bind most kinds of functions that rely on these types.
+One major limitation of the above is that every data conversion implicitly
-One major caveat are functions that take Eigen matrices *by reference* and modify
+involves a copy, which can be both expensive (for large matrices) and disallows
-them somehow, in which case the information won't be propagated to the caller.
+binding functions that change their (Matrix) arguments.  Pybind11 allows you to
+work around this by using Eigen's ``Eigen::Ref<MatrixType>`` class much as you
+would when writing a function taking a generic type in Eigen itself (subject to
+some limitations discussed below).
+
+When calling a bound function accepting a ``Eigen::Ref<const MatrixType>``
+type, pybind11 will attempt to avoid copying by using an ``Eigen::Map`` object
+that maps into the source ``numpy.ndarray`` data: this requires both that the
+data types are the same (e.g. ``dtype='float64'`` and ``MatrixType::Scalar`` is
+``double``); and that the storage is layout compatible.  The latter limitation
+is discussed in detail in the section below, and requires careful
+consideration: by default, numpy matrices and eigen matrices are *not* storage
+compatible.
+
+If the numpy matrix cannot be used as is (either because its types differ, e.g.
+passing an array of integers to an Eigen paramater requiring doubles, or
+because the storage is incompatible), pybind11 makes a temporary copy and
+passes the copy instead.
+
+When a bound function parameter is instead ``Eigen::Ref<MatrixType>`` (note the
+lack of ``const``), pybind11 will only allow the function to be called if it
+can be mapped *and* if the numpy array is writeable (that is
+``a.flags.writeable`` is true).  Any access (including modification) made to
+the passed variable will be transparently carried out directly on the
+``numpy.ndarray``.
+
+This means you can can write code such as the following and have it work as
+expected:
 
 .. code-block:: cpp
 
-    /* The Python bindings of these functions won't replicate
+    void scale_by_2(Eigen::Ref<Eigen::VectorXd> m) {
-       the intended effect of modifying the function arguments */
-    void scale_by_2(Eigen::Vector3f &v) {
-        v *= 2;
-    }
-    void scale_by_2(Eigen::Ref<Eigen::MatrixXd> &v) {
         v *= 2;
     }
 
-To see why this is, refer to the section on :ref:`opaque` (although that
+Note, however, that you will likely run into limitations due to numpy and
-section specifically covers STL data types, the underlying issue is the same).
+Eigen's difference default storage order for data; see the below section on
-The :ref:`numpy` sections discuss an efficient alternative for exposing the
+:ref:`storage_orders` for details on how to bind code that won't run into such
-underlying native Eigen types as opaque objects in a way that still integrates
+limitations.
-with NumPy and SciPy.
+
+.. note::
+
+    Passing by reference is not supported for sparse types.
+
+Returning values to Python
+==========================
+
+When returning an ordinary dense Eigen matrix type to numpy (e.g.
+``Eigen::MatrixXd`` or ``Eigen::RowVectorXf``) pybind11 keeps the matrix and
+returns a numpy array that directly references the Eigen matrix: no copy of the
+data is performed.  The numpy array will have ``array.flags.owndata`` set to
+``False`` to indicate that it does not own the data, and the lifetime of the
+stored Eigen matrix will be tied to the returned ``array``.
+
+If you bind a function with a non-reference, ``const`` return type (e.g.
+``const Eigen::MatrixXd``), the same thing happens except that pybind11 also
+sets the numpy array's ``writeable`` flag to false.
+
+If you return an lvalue reference or pointer, the usual pybind11 rules apply,
+as dictated by the binding function's return value policy (see the
+documentation on :ref:`return_value_policies` for full details).  That means,
+without an explicit return value policy, lvalue references will be copied and
+pointers will be managed by pybind11.  In order to avoid copying, you should
+explictly specify an appropriate return value policy, as in the following
+example:
+
+.. code-block:: cpp
+
+    class MyClass {
+        Eigen::MatrixXd big_mat = Eigen::MatrixXd::Zero(10000, 10000);
+    public:
+        Eigen::MatrixXd &getMatrix() { return big_mat; }
+        const Eigen::MatrixXd &viewMatrix() { return big_mat; }
+    };
+
+    // Later, in binding code:
+    py::class_<MyClass>(m, "MyClass")
+        .def(py::init<>())
+        .def("copy_matrix", &MyClass::getMatrix) // Makes a copy!
+        .def("get_matrix", &MyClass::getMatrix, py::return_value_policy::reference_internal)
+        .def("view_matrix", &MyClass::viewMatrix, py::return_value_policy::reference_internal)
+        ;
+
+.. code-block:: python
+
+    a = MyClass()
+    m = a.get_matrix()   # flags.writeable = True,  flags.owndata = False
+    v = a.view_matrix()  # flags.writeable = False, flags.owndata = False
+    c = a.copy_matrix()  # flags.writeable = True,  flags.owndata = True
+    # m[5,6] and v[5,6] refer to the same element, c[5,6] does not.
+
+Note in this example that ``py::return_value_policy::reference_internal`` is
+used to tie the life of the MyClass object to the life of the returned arrays.
+
+You may also return an ``Eigen::Ref``, ``Eigen::Map`` or other map-like Eigen
+object (for example, the return value of ``matrix.block()`` and related
+methods) that map into a dense Eigen type.  When doing so, the default
+behaviour of pybind11 is to simply reference the returned data: you must take
+care to ensure that this data remains valid!  You may ask pybind11 to
+explicitly *copy* such a return value by using the
+``py::return_value_policy::copy`` policy when binding the function.  You may
+also use ``py::return_value_policy::reference_internal`` or a
+``py::keep_alive`` to ensure the data stays valid as long as the returned numpy
+array does.
+
+When returning such a reference of map, pybind11 additionally respects the
+readonly-status of the returned value, marking the numpy array as non-writeable
+if the reference or map was itself read-only.
+
+.. note::
+
+    Sparse types are always copied when returned.
+
+.. _storage_orders:
+
+Storage orders
+==============
+
+Passing arguments via ``Eigen::Ref`` has some limitations that you must be
+aware of in order to effectively pass matrices by reference.  First and
+foremost is that the default ``Eigen::Ref<MatrixType>`` class requires
+contiguous storage along columns (for column-major types, the default in Eigen)
+or rows if ``MatrixType`` is specifically an ``Eigen::RowMajor`` storage type.
+The former, Eigen's default, is incompatible with ``numpy``'s default row-major
+storage, and so you will not be able to pass numpy arrays to Eigen by reference
+without making one of two changes.
+
+(Note that this does not apply to vectors (or column or row matrices): for such
+types the "row-major" and "column-major" distinction is meaningless).
+
+The first approach is to change the use of ``Eigen::Ref<MatrixType>`` to the
+more general ``Eigen::Ref<MatrixType, 0, Eigen::Stride<Eigen::Dynamic,
+Eigen::Dynamic>>`` (or similar type with a fully dynamic stride type in the
+third template argument).  Since this is a rather cumbersome type, pybind11
+provides a ``py::EigenDRef<MatrixType>`` type alias for your convenience (along
+with EigenDMap for the equivalent Map, and EigenDStride for just the stride
+type).
+
+This type allows Eigen to map into any arbitrary storage order.  This is not
+the default in Eigen for performance reasons: contiguous storage allows
+vectorization that cannot be done when storage is not known to be contiguous at
+compile time.  The default ``Eigen::Ref`` stride type allows non-contiguous
+storage along the outer dimension (that is, the rows of a column-major matrix
+or columns of a row-major matrix), but not along the inner dimension.
+
+This type, however, has the added benefit of also being able to map numpy array
+slices.  For example, the following (contrived) example uses Eigen with a numpy
+slice to multiply by 2 all coefficients that are both on even rows (0, 2, 4,
+...) and in columns 2, 5, or 8:
+
+.. code-block:: cpp
+
+    m.def("scale", [](py::EigenDRef<Eigen::MatrixXd> m, double c) { m *= c; });
+
+.. code-block:: python
+
+    # a = np.array(...)
+    scale_by_2(myarray[0::2, 2:9:3])
+
+The second approach to avoid copying is more intrusive: rearranging the
+underlying data types to not run into the non-contiguous storage problem in the
+first place.  In particular, that means using matrices with ``Eigen::RowMajor``
+storage, where appropriate, such as:
+
+.. code-block:: cpp
+
+    using RowMatrixXd = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+    // Use RowMatrixXd instead of MatrixXd
+
+Now bound functions accepting ``Eigen::Ref<RowMatrixXd>`` arguments will be
+callable with numpy's (default) arrays without involving a copying.
+
+You can, alternatively, change the storage order that numpy arrays use by
+adding the ``order='F'`` option when creating an array:
+
+.. code-block:: python
+
+    myarray = np.array(source, order='F')
+
+Such an object will be passable to a bound function accepting an
+``Eigen::Ref<MatrixXd>`` (or similar column-major Eigen type).
+
+One major caveat with this approach, however, is that it is not entirely as
+easy as simply flipping all Eigen or numpy usage from one to the other: some
+operations may alter the storage order of a numpy array.  For example, ``a2 =
+array.transpose()`` results in ``a2`` being a view of ``array`` that references
+the same data, but in the opposite storage order!
+
+While this approach allows fully optimized vectorized calculations in Eigen, it
+cannot be used with array slices, unlike the first approach.
+
+When *returning* a matrix to Python (either a regular matrix, a reference via
+``Eigen::Ref<>``, or a map/block into a matrix), no special storage
+consideration is required: the created numpy array will have the required
+stride that allows numpy to properly interpret the array, whatever its storage
+order.
+
+Failing rather than copying
+===========================
+
+The default behaviour when binding ``Eigen::Ref<const MatrixType>`` eigen
+references is to copy matrix values when passed a numpy array that does not
+conform to the element type of ``MatrixType`` or does not have a compatible
+stride layout.  If you want to explicitly avoid copying in such a case, you
+should bind arguments using the ``py::arg().noconvert()`` annotation (as
+described in the :ref:`nonconverting_arguments` documentation).
+
+The following example shows an example of arguments that don't allow data
+copying to take place:
+
+.. code-block:: cpp
+
+    // The method and function to be bound:
+    class MyClass {
+        // ...
+        double some_method(const Eigen::Ref<const MatrixXd> &matrix) { /* ... */ }
+    };
+    float some_function(const Eigen::Ref<const MatrixXf> &big,
+                        const Eigen::Ref<const MatrixXf> &small) {
+        // ...
+    }
+
+    // The associated binding code:
+    using namespace pybind11::literals; // for "arg"_a
+    py::class_<MyClass>(m, "MyClass")
+        // ... other class definitions
+        .def("some_method", &MyClass::some_method, py::arg().nocopy());
+
+    m.def("some_function", &some_function,
+        "big"_a.nocopy(), // <- Don't allow copying for this arg
+        "small"_a         // <- This one can be copied if needed
+    );
+
+With the above binding code, attempting to call the the ``some_method(m)``
+method on a ``MyClass`` object, or attempting to call ``some_function(m, m2)``
+will raise a ``RuntimeError`` rather than making a temporary copy of the array.
+It will, however, allow the ``m2`` argument to be copied into a temporary if
+necessary.
+
+Note that explicitly specifying ``.noconvert()`` is not required for *mutable*
+Eigen references (e.g. ``Eigen::Ref<MatrixXd>`` without ``const`` on the
+``MatrixXd``): mutable references will never be called with a temporary copy.
+
+Vectors versus column/row matrices
+==================================
+
+Eigen and numpy have fundamentally different notions of a vector.  In Eigen, a
+vector is simply a matrix with the number of columns or rows set to 1 at
+compile time (for a column vector or row vector, respectively).  Numpy, in
+contast, has comparable 2-dimensional 1xN and Nx1 arrays, but *also* has
+1-dimensional arrays of size N.
+
+When passing a 2-dimensional 1xN or Nx1 array to Eigen, the Eigen type must
+have matching dimensions: That is, you cannot pass a 2-dimensional Nx1 numpy
+array to an Eigen value expecting a row vector, or a 1xN numpy array as a
+column vector argument.
+
+On the other hand, pybind11 allows you to pass 1-dimensional arrays of length N
+as Eigen parameters.  If the Eigen type can hold a column vector of length N it
+will be passed as such a column vector.  If not, but the Eigen type constraints
+will accept a row vector, it will be passed as a row vector.  (The column
+vector takes precendence when both are supported, for example, when passing a
+1D numpy array to a MatrixXd argument).  Note that the type need not be
+expicitly a vector: it is permitted to pass a 1D numpy array of size 5 to an
+Eigen ``Matrix<double, Dynamic, 5>``: you would end up with a 1x5 Eigen matrix.
+Passing the same to an ``Eigen::MatrixXd`` would result in a 5x1 Eigen matrix.
+
+When returning an eigen vector to numpy, the conversion is ambiguous: a row
+vector of length 4 could be returned as either a 1D array of length 4, or as a
+2D array of size 1x4.  When encoutering such a situation, pybind11 compromises
+by considering the returned Eigen type: if it is a compile-time vector--that
+is, the type has either the number of rows or columns set to 1 at compile
+time--pybind11 converts to a 1D numpy array when returning the value.  For
+instances that are a vector only at run-time (e.g. ``MatrixXd``,
+``Matrix<float, Dynamic, 4>``), pybind11 returns the vector as a 2D array to
+numpy.  If this isn't want you want, you can use ``array.reshape(...)`` to get
+a view of the same data in the desired dimensions.
 
 .. seealso::
 

docs/advanced/functions.rst

@@ -6,6 +6,8 @@ with the basics of binding functions and classes, as explained in :doc:`/basics`
 and :doc:`/classes`. The following guide is applicable to both free and member
 functions, i.e. *methods* in Python.
 
+.. _return_value_policies:
+
 Return value policies
 =====================
 
@@ -320,6 +322,8 @@ like so:
     py::class_<MyClass>("MyClass")
         .def("myFunction", py::arg("arg") = (SomeType *) nullptr);
 
+.. _nonconverting_arguments:
+
 Non-converting arguments
 ========================
 

include/pybind11/eigen.h

@@ -30,156 +30,487 @@
 #  pragma warning(disable: 4127) // warning C4127: Conditional expression is constant
 #endif
 
+// Eigen prior to 3.2.7 doesn't have proper move constructors--but worse, some classes get implicit
+// move constructors that break things.  We could detect this an explicitly copy, but an extra copy
+// of matrices seems highly undesirable.
+static_assert(EIGEN_VERSION_AT_LEAST(3,2,7), "Eigen support in pybind11 requires Eigen >= 3.2.7");
+
 NAMESPACE_BEGIN(pybind11)
+
+// Provide a convenience alias for easier pass-by-ref usage with fully dynamic strides:
+using EigenDStride = Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>;
+template <typename MatrixType> using EigenDRef = Eigen::Ref<MatrixType, 0, EigenDStride>;
+template <typename MatrixType> using EigenDMap = Eigen::Map<MatrixType, 0, EigenDStride>;
+
 NAMESPACE_BEGIN(detail)
 
-template <typename T> using is_eigen_dense = is_template_base_of<Eigen::DenseBase, T>;
+#if EIGEN_VERSION_AT_LEAST(3,3,0)
-template <typename T> using is_eigen_sparse = is_template_base_of<Eigen::SparseMatrixBase, T>;
+using EigenIndex = Eigen::Index;
-template <typename T> using is_eigen_ref = is_template_base_of<Eigen::RefBase, T>;
+#else
+using EigenIndex = EIGEN_DEFAULT_DENSE_INDEX_TYPE;
+#endif
 
+// Matches Eigen::Map, Eigen::Ref, blocks, etc:
+template <typename T> using is_eigen_dense_map = all_of<is_template_base_of<Eigen::DenseBase, T>, std::is_base_of<Eigen::MapBase<T, Eigen::ReadOnlyAccessors>, T>>;
+template <typename T> using is_eigen_mutable_map = std::is_base_of<Eigen::MapBase<T, Eigen::WriteAccessors>, T>;
+template <typename T> using is_eigen_dense_plain = all_of<negation<is_eigen_dense_map<T>>, is_template_base_of<Eigen::PlainObjectBase, T>>;
+template <typename T> using is_eigen_sparse = is_template_base_of<Eigen::SparseMatrixBase, T>;
 // Test for objects inheriting from EigenBase<Derived> that aren't captured by the above.  This
 // basically covers anything that can be assigned to a dense matrix but that don't have a typical
 // matrix data layout that can be copied from their .data().  For example, DiagonalMatrix and
 // SelfAdjointView fall into this category.
-template <typename T> using is_eigen_base = all_of<
+template <typename T> using is_eigen_other = all_of<
     is_template_base_of<Eigen::EigenBase, T>,
-    negation<is_eigen_dense<T>>,
+    negation<any_of<is_eigen_dense_map<T>, is_eigen_dense_plain<T>, is_eigen_sparse<T>>>
-    negation<is_eigen_sparse<T>>
 >;
 
+// Captures numpy/eigen conformability status (returned by EigenProps::conformable()):
+template <bool EigenRowMajor> struct EigenConformable {
+    bool conformable = false;
+    EigenIndex rows = 0, cols = 0;
+    EigenDStride stride{0, 0};
+
+    EigenConformable(bool fits = false) : conformable{fits} {}
+    // Matrix type:
+    EigenConformable(EigenIndex r, EigenIndex c,
+            EigenIndex rstride, EigenIndex cstride) :
+        conformable{true}, rows{r}, cols{c},
+        stride(EigenRowMajor ? rstride : cstride /* outer stride */,
+               EigenRowMajor ? cstride : rstride /* inner stride */)
+        {}
+    // Vector type:
+    EigenConformable(EigenIndex r, EigenIndex c, EigenIndex stride) : EigenConformable(r, c, r == 1 ? c*stride : stride, c == 1 ? r : r*stride) {}
+    template <typename props> bool stride_compatible() const {
+        return
+            (props::inner_stride == Eigen::Dynamic || props::inner_stride == stride.inner()) &&
+            (props::outer_stride == Eigen::Dynamic || props::outer_stride == stride.outer());
+    }
+    operator bool() const { return conformable; }
+};
+
+template <typename Type> struct eigen_extract_stride { using type = Type; };
+template <typename PlainObjectType, int MapOptions, typename StrideType>
+struct eigen_extract_stride<Eigen::Map<PlainObjectType, MapOptions, StrideType>> { using type = StrideType; };
+template <typename PlainObjectType, int Options, typename StrideType>
+struct eigen_extract_stride<Eigen::Ref<PlainObjectType, Options, StrideType>> { using type = StrideType; };
+
+// Helper struct for extracting information from an Eigen type
+template <typename Type_> struct EigenProps {
+    using Type = Type_;
+    using Scalar = typename Type::Scalar;
+    using StrideType = typename eigen_extract_stride<Type>::type;
+    static constexpr EigenIndex
+        rows = Type::RowsAtCompileTime,
+        cols = Type::ColsAtCompileTime,
+        size = Type::SizeAtCompileTime;
+    static constexpr bool
+        row_major = Type::IsRowMajor,
+        vector = Type::IsVectorAtCompileTime, // At least one dimension has fixed size 1
+        fixed_rows = rows != Eigen::Dynamic,
+        fixed_cols = cols != Eigen::Dynamic,
+        fixed = size != Eigen::Dynamic, // Fully-fixed size
+        dynamic = !fixed_rows && !fixed_cols; // Fully-dynamic size
+
+    template <EigenIndex i, EigenIndex ifzero> using if_zero = std::integral_constant<EigenIndex, i == 0 ? ifzero : i>;
+    static constexpr EigenIndex inner_stride = if_zero<StrideType::InnerStrideAtCompileTime, 1>::value,
+                                outer_stride = if_zero<StrideType::OuterStrideAtCompileTime,
+                                                       vector ? size : row_major ? cols : rows>::value;
+    static constexpr bool dynamic_stride = inner_stride == Eigen::Dynamic && outer_stride == Eigen::Dynamic;
+    static constexpr bool requires_row_major = !dynamic_stride && !vector && (row_major ? inner_stride : outer_stride) == 1;
+    static constexpr bool requires_col_major = !dynamic_stride && !vector && (row_major ? outer_stride : inner_stride) == 1;
+
+    // Takes an input array and determines whether we can make it fit into the Eigen type.  If
+    // the array is a vector, we attempt to fit it into either an Eigen 1xN or Nx1 vector
+    // (preferring the latter if it will fit in either, i.e. for a fully dynamic matrix type).
+    static EigenConformable<row_major> conformable(const array &a) {
+        const auto dims = a.ndim();
+        if (dims < 1 || dims > 2)
+            return false;
+
+        if (dims == 2) { // Matrix type: require exact match (or dynamic)
+
+            EigenIndex
+                np_rows = a.shape(0),
+                np_cols = a.shape(1),
+                np_rstride = a.strides(0) / sizeof(Scalar),
+                np_cstride = a.strides(1) / sizeof(Scalar);
+            if ((fixed_rows && np_rows != rows) || (fixed_cols && np_cols != cols))
+                return false;
+
+            return {np_rows, np_cols, np_rstride, np_cstride};
+        }
+
+        // Otherwise we're storing an n-vector.  Only one of the strides will be used, but whichever
+        // is used, we want the (single) numpy stride value.
+        const EigenIndex n = a.shape(0),
+              stride = a.strides(0) / sizeof(Scalar);
+
+        if (vector) { // Eigen type is a compile-time vector
+            if (fixed && size != n)
+                return false; // Vector size mismatch
+            return {rows == 1 ? 1 : n, cols == 1 ? 1 : n, stride};
+        }
+        else if (fixed) {
+            // The type has a fixed size, but is not a vector: abort
+            return false;
+        }
+        else if (fixed_cols) {
+            // Since this isn't a vector, cols must be != 1.  We allow this only if it exactly
+            // equals the number of elements (rows is Dynamic, and so 1 row is allowed).
+            if (cols != n) return false;
+            return {1, n, stride};
+        }
+        else {
+            // Otherwise it's either fully dynamic, or column dynamic; both become a column vector
+            if (fixed_rows && rows != n) return false;
+            return {n, 1, stride};
+        }
+    }
+
+    static PYBIND11_DESCR descriptor() {
+        constexpr bool show_writeable = is_eigen_dense_map<Type>::value && is_eigen_mutable_map<Type>::value;
+        constexpr bool show_order = is_eigen_dense_map<Type>::value;
+        constexpr bool show_c_contiguous = show_order && requires_row_major;
+        constexpr bool show_f_contiguous = !show_c_contiguous && show_order && requires_col_major;
+
+    return _("numpy.ndarray[") + npy_format_descriptor<Scalar>::name() +
+        _("[")  + _<fixed_rows>(_<(size_t) rows>(), _("m")) +
+        _(", ") + _<fixed_cols>(_<(size_t) cols>(), _("n")) +
+        _("]") +
+        // For a reference type (e.g. Ref<MatrixXd>) we have other constraints that might need to be
+        // satisfied: writeable=True (for a mutable reference), and, depending on the map's stride
+        // options, possibly f_contiguous or c_contiguous.  We include them in the descriptor output
+        // to provide some hint as to why a TypeError is occurring (otherwise it can be confusing to
+        // see that a function accepts a 'numpy.ndarray[float64[3,2]]' and an error message that you
+        // *gave* a numpy.ndarray of the right type and dimensions.
+        _<show_writeable>(", flags.writeable", "") +
+        _<show_c_contiguous>(", flags.c_contiguous", "") +
+        _<show_f_contiguous>(", flags.f_contiguous", "") +
+        _("]");
+    }
+};
+
+// Casts an Eigen type to numpy array.  If given a base, the numpy array references the src data,
+// otherwise it'll make a copy.  writeable lets you turn off the writeable flag for the array.
+template <typename props> handle eigen_array_cast(typename props::Type const &src, handle base = handle(), bool writeable = true) {
+    constexpr size_t elem_size = sizeof(typename props::Scalar);
+    std::vector<size_t> shape, strides;
+    if (props::vector) {
+        shape.push_back(src.size());
+        strides.push_back(elem_size * src.innerStride());
+    }
+    else {
+        shape.push_back(src.rows());
+        shape.push_back(src.cols());
+        strides.push_back(elem_size * src.rowStride());
+        strides.push_back(elem_size * src.colStride());
+    }
+    array a(std::move(shape), std::move(strides), src.data(), base);
+    if (!writeable)
+        array_proxy(a.ptr())->flags &= ~detail::npy_api::NPY_ARRAY_WRITEABLE_;
+
+    return a.release();
+}
+
+// Takes an lvalue ref to some Eigen type and a (python) base object, creating a numpy array that
+// reference the Eigen object's data with `base` as the python-registered base class (if omitted,
+// the base will be set to None, and lifetime management is up to the caller).  The numpy array is
+// non-writeable if the given type is const.
+template <typename props, typename Type>
+handle eigen_ref_array(Type &src, handle parent = none()) {
+    // none here is to get past array's should-we-copy detection, which currently always
+    // copies when there is no base.  Setting the base to None should be harmless.
+    return eigen_array_cast<props>(src, parent, !std::is_const<Type>::value);
+}
+
+// Takes a pointer to some dense, plain Eigen type, builds a capsule around it, then returns a numpy
+// array that references the encapsulated data with a python-side reference to the capsule to tie
+// its destruction to that of any dependent python objects.  Const-ness is determined by whether or
+// not the Type of the pointer given is const.
+template <typename props, typename Type, typename = enable_if_t<is_eigen_dense_plain<Type>::value>>
+handle eigen_encapsulate(Type *src) {
+    capsule base(src, [](PyObject *o) { delete reinterpret_steal<capsule>(o).operator Type*(); });
+    return eigen_ref_array<props>(*src, base);
+}
+
+// Type caster for regular, dense matrix types (e.g. MatrixXd), but not maps/refs/etc. of dense
+// types.
 template<typename Type>
-struct type_caster<Type, enable_if_t<is_eigen_dense<Type>::value && !is_eigen_ref<Type>::value>> {
+struct type_caster<Type, enable_if_t<is_eigen_dense_plain<Type>::value>> {
-    typedef typename Type::Scalar Scalar;
+    using Scalar = typename Type::Scalar;
-    static constexpr bool rowMajor = Type::IsRowMajor;
+    using props = EigenProps<Type>;
-    static constexpr bool isVector = Type::IsVectorAtCompileTime;
-    static constexpr auto nRows = Type::RowsAtCompileTime, nCols = Type::ColsAtCompileTime;
 
     bool load(handle src, bool) {
         auto buf = array_t<Scalar>::ensure(src);
         if (!buf)
             return false;
 
-        using namespace Eigen;
+        auto dims = buf.ndim();
-
+        if (dims < 1 || dims > 2)
-        using Strides = Stride<Dynamic, Dynamic>;
-        using AMap = Map<Type, 0, Strides>;
-
-        if (buf.ndim() == 1) { // A one-dimensional array
-            Index n_elts = buf.shape(0);
-            size_t str = buf.strides(0) / sizeof(Scalar);
-            Strides stride(str, str); // Whether we map to inner or outer is irrelevant
-            if (isVector) {
-                if (Type::SizeAtCompileTime != Dynamic && Type::SizeAtCompileTime != n_elts)
-                    return false; // Vector size mismatch
-                value = AMap(buf.mutable_data(), nRows == 1 ? 1 : n_elts, nCols == 1 ? 1 : n_elts, stride);
-            }
-            else if (Type::SizeAtCompileTime != Dynamic) {
-                // The type has a fixed size, but is not a vector: abort
-                return false;
-            }
-            else if (nRows == Dynamic && nCols == Dynamic) {
-                // Fully dynamic size.  numpy doesn't distinguish between a row vector and column
-                // vector, so we'll (arbitrarily) choose a column vector.
-                value = AMap(buf.mutable_data(), n_elts, 1, stride);
-            }
-            else if (nRows != Dynamic) {
-                // Since this isn't a vector, nRows must be != 1.  We allow this only if it exactly
-                // equals the number of elements (nCols is Dynamic, and so 1 is allowed).
-                if (nRows != n_elts) return false;
-                value = AMap(buf.mutable_data(), n_elts, 1, stride);
-            }
-            else { // nCols != Dynamic; same as above, but for fixed columns
-                if (nCols != n_elts) return false;
-                value = AMap(buf.mutable_data(), 1, n_elts, stride);
-            }
-        } else if (buf.ndim() == 2) {
-
-            if ((Type::RowsAtCompileTime != Dynamic && buf.shape(0) != (size_t) Type::RowsAtCompileTime) ||
-                (Type::ColsAtCompileTime != Dynamic && buf.shape(1) != (size_t) Type::ColsAtCompileTime))
             return false;
 
-            value = AMap(
+        auto fits = props::conformable(buf);
-                buf.mutable_data(), buf.shape(0), buf.shape(1),
+        if (!fits)
-                Strides(buf.strides(rowMajor ? 0 : 1) / sizeof(Scalar),
+            return false; // Non-comformable vector/matrix types
-                        buf.strides(rowMajor ? 1 : 0) / sizeof(Scalar))
+
-            );
+        value = Eigen::Map<const Type, 0, EigenDStride>(buf.data(), fits.rows, fits.cols, fits.stride);
-        } else {
+
-            return false;
-        }
         return true;
     }
 
-    static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) {
+private:
-        if (isVector) {
-            return array(
-                { (size_t) src.size() },                                      // shape
-                { sizeof(Scalar) * static_cast<size_t>(src.innerStride()) },  // strides
-                src.data()                                                    // data
-            ).release();
-        } else {
-            return array(
-                { (size_t) src.rows(),                                        // shape
-                  (size_t) src.cols() },
-                { sizeof(Scalar) * static_cast<size_t>(src.rowStride()),      // strides
-                  sizeof(Scalar) * static_cast<size_t>(src.colStride()) },
-                src.data()                                                    // data
-            ).release();
-        }
-    }
 
-    PYBIND11_TYPE_CASTER(Type, _("numpy.ndarray[") + npy_format_descriptor<Scalar>::name() +
+    // Cast implementation
-            _("[") + rows() + _(", ") + cols() + _("]]"));
+    template <typename CType>
-
+    static handle cast_impl(CType *src, return_value_policy policy, handle parent) {
-protected:
+        switch (policy) {
-    template <typename T = Type, enable_if_t<T::RowsAtCompileTime == Eigen::Dynamic, int> = 0>
+            case return_value_policy::take_ownership:
-    static PYBIND11_DESCR rows() { return _("m"); }
+            case return_value_policy::automatic:
-    template <typename T = Type, enable_if_t<T::RowsAtCompileTime != Eigen::Dynamic, int> = 0>
+                return eigen_encapsulate<props>(src);
-    static PYBIND11_DESCR rows() { return _<T::RowsAtCompileTime>(); }
+            case return_value_policy::move:
-    template <typename T = Type, enable_if_t<T::ColsAtCompileTime == Eigen::Dynamic, int> = 0>
+                return eigen_encapsulate<props>(new CType(std::move(*src)));
-    static PYBIND11_DESCR cols() { return _("n"); }
+            case return_value_policy::copy:
-    template <typename T = Type, enable_if_t<T::ColsAtCompileTime != Eigen::Dynamic, int> = 0>
+                return eigen_array_cast<props>(*src);
-    static PYBIND11_DESCR cols() { return _<T::ColsAtCompileTime>(); }
+            case return_value_policy::reference:
+            case return_value_policy::automatic_reference:
+                return eigen_ref_array<props>(*src);
+            case return_value_policy::reference_internal:
+                return eigen_ref_array<props>(*src, parent);
+            default:
+                throw cast_error("unhandled return_value_policy: should not happen!");
         };
+    }
 
-// Eigen::Ref<Derived> satisfies is_eigen_dense, but isn't constructable, so it needs a special
-// type_caster to handle argument copying/forwarding.
-template <typename CVDerived, int Options, typename StrideType>
-struct type_caster<Eigen::Ref<CVDerived, Options, StrideType>> {
-protected:
-    using Type = Eigen::Ref<CVDerived, Options, StrideType>;
-    using Derived = typename std::remove_const<CVDerived>::type;
-    using DerivedCaster = make_caster<Derived>;
-    DerivedCaster derived_caster;
-    std::unique_ptr<Type> value;
 public:
-    bool load(handle src, bool convert) { if (derived_caster.load(src, convert)) { value.reset(new Type(derived_caster.operator Derived&())); return true; } return false; }
-    static handle cast(const Type &src, return_value_policy policy, handle parent) { return DerivedCaster::cast(src, policy, parent); }
-    static handle cast(const Type *src, return_value_policy policy, handle parent) { return DerivedCaster::cast(*src, policy, parent); }
 
-    static PYBIND11_DESCR name() { return DerivedCaster::name(); }
+    // Normal returned non-reference, non-const value:
+    static handle cast(Type &&src, return_value_policy /* policy */, handle parent) {
+        return cast_impl(&src, return_value_policy::move, parent);
+    }
+    // If you return a non-reference const, we mark the numpy array readonly:
+    static handle cast(const Type &&src, return_value_policy /* policy */, handle parent) {
+        return cast_impl(&src, return_value_policy::move, parent);
+    }
+    // lvalue reference return; default (automatic) becomes copy
+    static handle cast(Type &src, return_value_policy policy, handle parent) {
+        if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference)
+            policy = return_value_policy::copy;
+        return cast_impl(&src, policy, parent);
+    }
+    // const lvalue reference return; default (automatic) becomes copy
+    static handle cast(const Type &src, return_value_policy policy, handle parent) {
+        if (policy == return_value_policy::automatic || policy == return_value_policy::automatic_reference)
+            policy = return_value_policy::copy;
+        return cast(&src, policy, parent);
+    }
+    // non-const pointer return
+    static handle cast(Type *src, return_value_policy policy, handle parent) {
+        return cast_impl(src, policy, parent);
+    }
+    // const pointer return
+    static handle cast(const Type *src, return_value_policy policy, handle parent) {
+        return cast_impl(src, policy, parent);
+    }
 
-    operator Type*() { return value.get(); }
+    static PYBIND11_DESCR name() { return type_descr(props::descriptor()); }
-    operator Type&() { if (!value) pybind11_fail("Eigen::Ref<...> value not loaded"); return *value; }
+
-    template <typename _T> using cast_op_type = pybind11::detail::cast_op_type<_T>;
+    operator Type*() { return &value; }
+    operator Type&() { return value; }
+    template <typename T> using cast_op_type = cast_op_type<T>;
+
+private:
+    Type value;
 };
 
-// type_caster for special matrix types (e.g. DiagonalMatrix): load() is not supported, but we can
+// Eigen Ref/Map classes have slightly different policy requirements, meaning we don't want to force
-// cast them into the python domain by first copying to a regular Eigen::Matrix, then casting that.
+// `move` when a Ref/Map rvalue is returned; we treat Ref<> sort of like a pointer (we care about
+// the underlying data, not the outer shell).
+template <typename Return>
+struct return_value_policy_override<Return, enable_if_t<is_eigen_dense_map<Return>::value>> {
+    static return_value_policy policy(return_value_policy p) { return p; }
+};
+
+// Base class for casting reference/map/block/etc. objects back to python.
+template <typename MapType> struct eigen_map_caster {
+private:
+    using props = EigenProps<MapType>;
+
+public:
+
+    // Directly referencing a ref/map's data is a bit dangerous (whatever the map/ref points to has
+    // to stay around), but we'll allow it under the assumption that you know what you're doing (and
+    // have an appropriate keep_alive in place).  We return a numpy array pointing directly at the
+    // ref's data (The numpy array ends up read-only if the ref was to a const matrix type.) Note
+    // that this means you need to ensure you don't destroy the object in some other way (e.g. with
+    // an appropriate keep_alive, or with a reference to a statically allocated matrix).
+    static handle cast(const MapType &src, return_value_policy policy, handle parent) {
+        switch (policy) {
+            case return_value_policy::copy:
+                return eigen_array_cast<props>(src);
+            case return_value_policy::reference_internal:
+                return eigen_array_cast<props>(src, parent, is_eigen_mutable_map<MapType>::value);
+            case return_value_policy::reference:
+            case return_value_policy::automatic:
+            case return_value_policy::automatic_reference:
+                return eigen_array_cast<props>(src, none(), is_eigen_mutable_map<MapType>::value);
+            default:
+                // move, take_ownership don't make any sense for a ref/map:
+                pybind11_fail("Invalid return_value_policy for Eigen Map/Ref/Block type");
+        }
+    }
+
+    static PYBIND11_DESCR name() { return props::descriptor(); }
+
+    // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return
+    // types but not bound arguments).  We still provide them (with an explicitly delete) so that
+    // you end up here if you try anyway.
+    bool load(handle, bool) = delete;
+    operator MapType() = delete;
+    template <typename> using cast_op_type = MapType;
+};
+
+// We can return any map-like object (but can only load Refs, specialized next):
+template <typename Type> struct type_caster<Type, enable_if_t<is_eigen_dense_map<Type>::value>>
+    : eigen_map_caster<Type> {};
+
+// Loader for Ref<...> arguments.  See the documentation for info on how to make this work without
+// copying (it requires some extra effort in many cases).
+template <typename PlainObjectType, typename StrideType>
+struct type_caster<
+    Eigen::Ref<PlainObjectType, 0, StrideType>,
+    enable_if_t<is_eigen_dense_map<Eigen::Ref<PlainObjectType, 0, StrideType>>::value>
+> : public eigen_map_caster<Eigen::Ref<PlainObjectType, 0, StrideType>> {
+private:
+    using Type = Eigen::Ref<PlainObjectType, 0, StrideType>;
+    using props = EigenProps<Type>;
+    using Scalar = typename props::Scalar;
+    using MapType = Eigen::Map<PlainObjectType, 0, StrideType>;
+    using Array = array_t<Scalar, array::forcecast |
+                ((props::row_major ? props::inner_stride : props::outer_stride) == 1 ? array::c_style :
+                 (props::row_major ? props::outer_stride : props::inner_stride) == 1 ? array::f_style : 0)>;
+    static constexpr bool need_writeable = is_eigen_mutable_map<Type>::value;
+    // Delay construction (these have no default constructor)
+    std::unique_ptr<MapType> map;
+    std::unique_ptr<Type> ref;
+    // Our array.  When possible, this is just a numpy array pointing to the source data, but
+    // sometimes we can't avoid copying (e.g. input is not a numpy array at all, has an incompatible
+    // layout, or is an array of a type that needs to be converted).  Using a numpy temporary
+    // (rather than an Eigen temporary) saves an extra copy when we need both type conversion and
+    // storage order conversion.  (Note that we refuse to use this temporary copy when loading an
+    // argument for a Ref<M> with M non-const, i.e. a read-write reference).
+    Array copy_or_ref;
+public:
+    bool load(handle src, bool convert) {
+        // First check whether what we have is already an array of the right type.  If not, we can't
+        // avoid a copy (because the copy is also going to do type conversion).
+        bool need_copy = !isinstance<Array>(src);
+
+        EigenConformable<props::row_major> fits;
+        if (!need_copy) {
+            // We don't need a converting copy, but we also need to check whether the strides are
+            // compatible with the Ref's stride requirements
+            Array aref = reinterpret_borrow<Array>(src);
+
+            if (aref && (!need_writeable || aref.writeable())) {
+                fits = props::conformable(aref);
+                if (!fits) return false; // Incompatible dimensions
+                if (!fits.template stride_compatible<props>())
+                    need_copy = true;
+                else
+                    copy_or_ref = std::move(aref);
+            }
+            else {
+                need_copy = true;
+            }
+        }
+
+        if (need_copy) {
+            // We need to copy: If we need a mutable reference, or we're not supposed to convert
+            // (either because we're in the no-convert overload pass, or because we're explicitly
+            // instructed not to copy (via `py::arg().noconvert()`) we have to fail loading.
+            if (!convert || need_writeable) return false;
+
+            Array copy = Array::ensure(src);
+            if (!copy) return false;
+            fits = props::conformable(copy);
+            if (!fits || !fits.template stride_compatible<props>())
+                return false;
+            copy_or_ref = copy;
+        }
+
+        ref.reset();
+        map.reset(new MapType(data(copy_or_ref), fits.rows, fits.cols, make_stride(fits.stride.outer(), fits.stride.inner())));
+        ref.reset(new Type(*map));
+
+        return true;
+    }
+
+    operator Type*() { return ref.get(); }
+    operator Type&() { return *ref; }
+    template <typename _T> using cast_op_type = pybind11::detail::cast_op_type<_T>;
+
+private:
+    template <typename T = Type, enable_if_t<is_eigen_mutable_map<T>::value, int> = 0>
+    Scalar *data(Array &a) { return a.mutable_data(); }
+
+    template <typename T = Type, enable_if_t<!is_eigen_mutable_map<T>::value, int> = 0>
+    const Scalar *data(Array &a) { return a.data(); }
+
+    // Attempt to figure out a constructor of `Stride` that will work.
+    // If both strides are fixed, use a default constructor:
+    template <typename S> using stride_ctor_default = bool_constant<
+        S::InnerStrideAtCompileTime != Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic &&
+        std::is_default_constructible<S>::value>;
+    // Otherwise, if there is a two-index constructor, assume it is (outer,inner) like
+    // Eigen::Stride, and use it:
+    template <typename S> using stride_ctor_dual = bool_constant<
+        !stride_ctor_default<S>::value && std::is_constructible<S, EigenIndex, EigenIndex>::value>;
+    // Otherwise, if there is a one-index constructor, and just one of the strides is dynamic, use
+    // it (passing whichever stride is dynamic).
+    template <typename S> using stride_ctor_outer = bool_constant<
+        !any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value &&
+        S::OuterStrideAtCompileTime == Eigen::Dynamic && S::InnerStrideAtCompileTime != Eigen::Dynamic &&
+        std::is_constructible<S, EigenIndex>::value>;
+    template <typename S> using stride_ctor_inner = bool_constant<
+        !any_of<stride_ctor_default<S>, stride_ctor_dual<S>>::value &&
+        S::InnerStrideAtCompileTime == Eigen::Dynamic && S::OuterStrideAtCompileTime != Eigen::Dynamic &&
+        std::is_constructible<S, EigenIndex>::value>;
+
+    template <typename S = StrideType, enable_if_t<stride_ctor_default<S>::value, int> = 0>
+    static S make_stride(EigenIndex, EigenIndex) { return S(); }
+    template <typename S = StrideType, enable_if_t<stride_ctor_dual<S>::value, int> = 0>
+    static S make_stride(EigenIndex outer, EigenIndex inner) { return S(outer, inner); }
+    template <typename S = StrideType, enable_if_t<stride_ctor_outer<S>::value, int> = 0>
+    static S make_stride(EigenIndex outer, EigenIndex) { return S(outer); }
+    template <typename S = StrideType, enable_if_t<stride_ctor_inner<S>::value, int> = 0>
+    static S make_stride(EigenIndex, EigenIndex inner) { return S(inner); }
+
+};
+
+// type_caster for special matrix types (e.g. DiagonalMatrix), which are EigenBase, but not
+// EigenDense (i.e. they don't have a data(), at least not with the usual matrix layout).
+// load() is not supported, but we can cast them into the python domain by first copying to a
+// regular Eigen::Matrix, then casting that.
 template <typename Type>
-struct type_caster<Type, enable_if_t<is_eigen_base<Type>::value && !is_eigen_ref<Type>::value>> {
+struct type_caster<Type, enable_if_t<is_eigen_other<Type>::value>> {
 protected:
-    using Matrix = Eigen::Matrix<typename Type::Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+    using Matrix = Eigen::Matrix<typename Type::Scalar, Type::RowsAtCompileTime, Type::ColsAtCompileTime>;
-    using MatrixCaster = make_caster<Matrix>;
+    using props = EigenProps<Matrix>;
 public:
-    [[noreturn]] bool load(handle, bool) { pybind11_fail("Unable to load() into specialized EigenBase object"); }
+    static handle cast(const Type &src, return_value_policy /* policy */, handle /* parent */) {
-    static handle cast(const Type &src, return_value_policy policy, handle parent) { return MatrixCaster::cast(Matrix(src), policy, parent); }
+        handle h = eigen_encapsulate<props>(new Matrix(src));
-    static handle cast(const Type *src, return_value_policy policy, handle parent) { return MatrixCaster::cast(Matrix(*src), policy, parent); }
+        return h;
+    }
+    static handle cast(const Type *src, return_value_policy policy, handle parent) { return cast(*src, policy, parent); }
 
-    static PYBIND11_DESCR name() { return MatrixCaster::name(); }
+    static PYBIND11_DESCR name() { return props::descriptor(); }
 
-    [[noreturn]] operator Type*() { pybind11_fail("Loading not supported for specialized EigenBase object"); }
+    // Explicitly delete these: support python -> C++ conversion on these (i.e. these can be return
-    [[noreturn]] operator Type&() { pybind11_fail("Loading not supported for specialized EigenBase object"); }
+    // types but not bound arguments).  We still provide them (with an explicitly delete) so that
-    template <typename _T> using cast_op_type = pybind11::detail::cast_op_type<_T>;
+    // you end up here if you try anyway.
+    bool load(handle, bool) = delete;
+    operator Type() = delete;
+    template <typename> using cast_op_type = Type;
 };
 
 template<typename Type>

tests/test_eigen.cpp

@@ -8,29 +8,47 @@
 */
 
 #include "pybind11_tests.h"
+#include "constructor_stats.h"
 #include <pybind11/eigen.h>
 #include <Eigen/Cholesky>
 
-Eigen::VectorXf double_col(const Eigen::VectorXf& x)
+using MatrixXdR = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
-{ return 2.0f * x; }
 
-Eigen::RowVectorXf double_row(const Eigen::RowVectorXf& x)
-{ return 2.0f * x; }
 
-Eigen::MatrixXf double_mat_cm(const Eigen::MatrixXf& x)
-{ return 2.0f * x; }
 
-// Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended
+// Sets/resets a testing reference matrix to have values of 10*r + c, where r and c are the
-Eigen::MatrixXd cholesky1(Eigen::Ref<Eigen::MatrixXd> &x) { return x.llt().matrixL(); }
+// (1-based) row/column number.
-Eigen::MatrixXd cholesky2(const Eigen::Ref<const Eigen::MatrixXd> &x) { return x.llt().matrixL(); }
+template <typename M> void reset_ref(M &x) {
-Eigen::MatrixXd cholesky3(const Eigen::Ref<Eigen::MatrixXd> &x) { return x.llt().matrixL(); }
+    for (int i = 0; i < x.rows(); i++) for (int j = 0; j < x.cols(); j++)
-Eigen::MatrixXd cholesky4(Eigen::Ref<const Eigen::MatrixXd> &x) { return x.llt().matrixL(); }
+        x(i, j) = 11 + 10*i + j;
-Eigen::MatrixXd cholesky5(Eigen::Ref<Eigen::MatrixXd> x) { return x.llt().matrixL(); }
+}
-Eigen::MatrixXd cholesky6(Eigen::Ref<const Eigen::MatrixXd> x) { return x.llt().matrixL(); }
 
-typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixXfRowMajor;
+// Returns a static, column-major matrix
-MatrixXfRowMajor double_mat_rm(const MatrixXfRowMajor& x)
+Eigen::MatrixXd &get_cm() {
-{ return 2.0f * x; }
+    static Eigen::MatrixXd *x;
+    if (!x) {
+        x = new Eigen::MatrixXd(3, 3);
+        reset_ref(*x);
+    }
+    return *x;
+}
+// Likewise, but row-major
+MatrixXdR &get_rm() {
+    static MatrixXdR *x;
+    if (!x) {
+        x = new MatrixXdR(3, 3);
+        reset_ref(*x);
+    }
+    return *x;
+}
+// Resets the values of the static matrices returned by get_cm()/get_rm()
+void reset_refs() {
+    reset_ref(get_cm());
+    reset_ref(get_rm());
+}
+
+// Returns element 2,1 from a matrix (used to test copy/nocopy)
+double get_elem(Eigen::Ref<const Eigen::MatrixXd> m) { return m(2, 1); };
 
 test_initializer eigen([](py::module &m) {
     typedef Eigen::Matrix<float, 5, 6, Eigen::RowMajor> FixedMatrixR;
@@ -46,21 +64,84 @@ test_initializer eigen([](py::module &m) {
 
     m.attr("have_eigen") = true;
 
-    // Non-symmetric matrix with zero elements
+    m.def("double_col", [](const Eigen::VectorXf &x) -> Eigen::VectorXf { return 2.0f * x; });
-    Eigen::MatrixXf mat(5, 6);
+    m.def("double_row", [](const Eigen::RowVectorXf &x) -> Eigen::RowVectorXf { return 2.0f * x; });
-    mat << 0, 3, 0, 0, 0, 11, 22, 0, 0, 0, 17, 11, 7, 5, 0, 1, 0, 11, 0,
+    m.def("double_threec", [](py::EigenDRef<Eigen::Vector3f> x) { x *= 2; });
-        0, 0, 0, 0, 11, 0, 0, 14, 0, 8, 11;
+    m.def("double_threer", [](py::EigenDRef<Eigen::RowVector3f> x) { x *= 2; });
+    m.def("double_mat_cm", [](Eigen::MatrixXf x) -> Eigen::MatrixXf { return 2.0f * x; });
+    m.def("double_mat_rm", [](DenseMatrixR x) -> DenseMatrixR { return 2.0f * x; });
 
-    m.def("double_col", &double_col);
+    // Different ways of passing via Eigen::Ref; the first and second are the Eigen-recommended
-    m.def("double_row", &double_row);
+    m.def("cholesky1", [](Eigen::Ref<MatrixXdR> x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
-    m.def("double_mat_cm", &double_mat_cm);
+    m.def("cholesky2", [](const Eigen::Ref<const MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
-    m.def("double_mat_rm", &double_mat_rm);
+    m.def("cholesky3", [](const Eigen::Ref<MatrixXdR> &x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
-    m.def("cholesky1", &cholesky1);
+    m.def("cholesky4", [](Eigen::Ref<const MatrixXdR> x) -> Eigen::MatrixXd { return x.llt().matrixL(); });
-    m.def("cholesky2", &cholesky2);
+
-    m.def("cholesky3", &cholesky3);
+    // Mutators: these add some value to the given element using Eigen, but Eigen should be mapping into
-    m.def("cholesky4", &cholesky4);
+    // the numpy array data and so the result should show up there.  There are three versions: one that
-    m.def("cholesky5", &cholesky5);
+    // works on a contiguous-row matrix (numpy's default), one for a contiguous-column matrix, and one
-    m.def("cholesky6", &cholesky6);
+    // for any matrix.
+    auto add_rm = [](Eigen::Ref<MatrixXdR> x, int r, int c, double v) { x(r,c) += v; };
+    auto add_cm = [](Eigen::Ref<Eigen::MatrixXd> x, int r, int c, double v) { x(r,c) += v; };
+
+    // Mutators (Eigen maps into numpy variables):
+    m.def("add_rm", add_rm); // Only takes row-contiguous
+    m.def("add_cm", add_cm); // Only takes column-contiguous
+    // Overloaded versions that will accept either row or column contiguous:
+    m.def("add1", add_rm);
+    m.def("add1", add_cm);
+    m.def("add2", add_cm);
+    m.def("add2", add_rm);
+    // This one accepts a matrix of any stride:
+    m.def("add_any", [](py::EigenDRef<Eigen::MatrixXd> x, int r, int c, double v) { x(r,c) += v; });
+
+    // Return mutable references (numpy maps into eigen varibles)
+    m.def("get_cm_ref", []() { return Eigen::Ref<Eigen::MatrixXd>(get_cm()); });
+    m.def("get_rm_ref", []() { return Eigen::Ref<MatrixXdR>(get_rm()); });
+    // The same references, but non-mutable (numpy maps into eigen variables, but is !writeable)
+    m.def("get_cm_const_ref", []() { return Eigen::Ref<const Eigen::MatrixXd>(get_cm()); });
+    m.def("get_rm_const_ref", []() { return Eigen::Ref<const MatrixXdR>(get_rm()); });
+    // Just the corners (via a Map instead of a Ref):
+    m.def("get_cm_corners", []() {
+        auto &x = get_cm();
+        return py::EigenDMap<Eigen::Matrix2d>(
+                x.data(),
+                py::EigenDStride(x.outerStride() * (x.rows() - 1), x.innerStride() * (x.cols() - 1)));
+    });
+    m.def("get_cm_corners_const", []() {
+        const auto &x = get_cm();
+        return py::EigenDMap<const Eigen::Matrix2d>(
+                x.data(),
+                py::EigenDStride(x.outerStride() * (x.rows() - 1), x.innerStride() * (x.cols() - 1)));
+    });
+
+    m.def("reset_refs", reset_refs); // Restores get_{cm,rm}_ref to original values
+
+    // Increments and returns ref to (same) matrix
+    m.def("incr_matrix", [](Eigen::Ref<Eigen::MatrixXd> m, double v) {
+        m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v);
+        return m;
+    }, py::return_value_policy::reference);
+
+    // Same, but accepts a matrix of any strides
+    m.def("incr_matrix_any", [](py::EigenDRef<Eigen::MatrixXd> m, double v) {
+        m += Eigen::MatrixXd::Constant(m.rows(), m.cols(), v);
+        return m;
+    }, py::return_value_policy::reference);
+
+    // Returns an eigen slice of even rows
+    m.def("even_rows", [](py::EigenDRef<Eigen::MatrixXd> m) {
+        return py::EigenDMap<Eigen::MatrixXd>(
+                m.data(), (m.rows() + 1) / 2, m.cols(),
+                py::EigenDStride(m.outerStride(), 2 * m.innerStride()));
+    }, py::return_value_policy::reference);
+
+    // Returns an eigen slice of even columns
+    m.def("even_cols", [](py::EigenDRef<Eigen::MatrixXd> m) {
+        return py::EigenDMap<Eigen::MatrixXd>(
+                m.data(), m.rows(), (m.cols() + 1) / 2,
+                py::EigenDStride(2 * m.outerStride(), m.innerStride()));
+    }, py::return_value_policy::reference);
 
     // Returns diagonals: a vector-like object with an inner stride != 1
     m.def("diagonal", [](const Eigen::Ref<const Eigen::MatrixXd> &x) { return x.diagonal(); });
@@ -72,6 +153,52 @@ test_initializer eigen([](py::module &m) {
         return x.block(start_row, start_col, block_rows, block_cols);
     });
 
+    // return value referencing/copying tests:
+    class ReturnTester {
+        Eigen::MatrixXd mat = create();
+    public:
+        ReturnTester() { print_created(this); }
+        ~ReturnTester() { print_destroyed(this); }
+        static Eigen::MatrixXd create() { return Eigen::MatrixXd::Ones(10, 10); }
+        static const Eigen::MatrixXd createConst() { return Eigen::MatrixXd::Ones(10, 10); }
+        Eigen::MatrixXd &get() { return mat; }
+        Eigen::MatrixXd *getPtr() { return &mat; }
+        const Eigen::MatrixXd &view() { return mat; }
+        const Eigen::MatrixXd *viewPtr() { return &mat; }
+        Eigen::Ref<Eigen::MatrixXd> ref() { return mat; }
+        Eigen::Ref<const Eigen::MatrixXd> refConst() { return mat; }
+        Eigen::Block<Eigen::MatrixXd> block(int r, int c, int nrow, int ncol) { return mat.block(r, c, nrow, ncol); }
+        Eigen::Block<const Eigen::MatrixXd> blockConst(int r, int c, int nrow, int ncol) const { return mat.block(r, c, nrow, ncol); }
+        py::EigenDMap<Eigen::Matrix2d> corners() { return py::EigenDMap<Eigen::Matrix2d>(mat.data(),
+                    py::EigenDStride(mat.outerStride() * (mat.outerSize()-1), mat.innerStride() * (mat.innerSize()-1))); }
+        py::EigenDMap<const Eigen::Matrix2d> cornersConst() const { return py::EigenDMap<const Eigen::Matrix2d>(mat.data(),
+                    py::EigenDStride(mat.outerStride() * (mat.outerSize()-1), mat.innerStride() * (mat.innerSize()-1))); }
+    };
+    using rvp = py::return_value_policy;
+    py::class_<ReturnTester>(m, "ReturnTester")
+        .def(py::init<>())
+        .def_static("create", &ReturnTester::create)
+        .def_static("create_const", &ReturnTester::createConst)
+        .def("get", &ReturnTester::get, rvp::reference_internal)
+        .def("get_ptr", &ReturnTester::getPtr, rvp::reference_internal)
+        .def("view", &ReturnTester::view, rvp::reference_internal)
+        .def("view_ptr", &ReturnTester::view, rvp::reference_internal)
+        .def("copy_get", &ReturnTester::get)   // Default rvp: copy
+        .def("copy_view", &ReturnTester::view) //         "
+        .def("ref", &ReturnTester::ref) // Default for Ref is to reference
+        .def("ref_const", &ReturnTester::refConst) // Likewise, but const
+        .def("ref_safe", &ReturnTester::ref, rvp::reference_internal)
+        .def("ref_const_safe", &ReturnTester::refConst, rvp::reference_internal)
+        .def("copy_ref", &ReturnTester::ref, rvp::copy)
+        .def("copy_ref_const", &ReturnTester::refConst, rvp::copy)
+        .def("block", &ReturnTester::block)
+        .def("block_safe", &ReturnTester::block, rvp::reference_internal)
+        .def("block_const", &ReturnTester::blockConst, rvp::reference_internal)
+        .def("copy_block", &ReturnTester::block, rvp::copy)
+        .def("corners", &ReturnTester::corners, rvp::reference_internal)
+        .def("corners_const", &ReturnTester::cornersConst, rvp::reference_internal)
+        ;
+
     // Returns a DiagonalMatrix with diagonal (1,2,3,...)
     m.def("incr_diag", [](int k) {
         Eigen::DiagonalMatrix<int, Eigen::Dynamic> m(k);
@@ -88,64 +215,49 @@ test_initializer eigen([](py::module &m) {
             return m.selfadjointView<Eigen::Upper>();
     });
 
-    m.def("fixed_r", [mat]() -> FixedMatrixR {
+    // Test matrix for various functions below.
-        return FixedMatrixR(mat);
+    Eigen::MatrixXf mat(5, 6);
-    });
+    mat << 0,  3,  0,  0,  0, 11,
+           22, 0,  0,  0, 17, 11,
+           7,  5,  0,  1,  0, 11,
+           0,  0,  0,  0,  0, 11,
+           0,  0, 14,  0,  8, 11;
 
-    m.def("fixed_c", [mat]() -> FixedMatrixC {
+    m.def("fixed_r", [mat]() -> FixedMatrixR { return FixedMatrixR(mat); });
-        return FixedMatrixC(mat);
+    m.def("fixed_r_const", [mat]() -> const FixedMatrixR { return FixedMatrixR(mat); });
-    });
+    m.def("fixed_c", [mat]() -> FixedMatrixC { return FixedMatrixC(mat); });
+    m.def("fixed_copy_r", [](const FixedMatrixR &m) -> FixedMatrixR { return m; });
+    m.def("fixed_copy_c", [](const FixedMatrixC &m) -> FixedMatrixC { return m; });
+    m.def("fixed_mutator_r", [](Eigen::Ref<FixedMatrixR>) {});
+    m.def("fixed_mutator_c", [](Eigen::Ref<FixedMatrixC>) {});
+    m.def("fixed_mutator_a", [](py::EigenDRef<FixedMatrixC>) {});
+    m.def("dense_r", [mat]() -> DenseMatrixR { return DenseMatrixR(mat); });
+    m.def("dense_c", [mat]() -> DenseMatrixC { return DenseMatrixC(mat); });
+    m.def("dense_copy_r", [](const DenseMatrixR &m) -> DenseMatrixR { return m; });
+    m.def("dense_copy_c", [](const DenseMatrixC &m) -> DenseMatrixC { return m; });
+    m.def("sparse_r", [mat]() -> SparseMatrixR { return Eigen::SparseView<Eigen::MatrixXf>(mat); });
+    m.def("sparse_c", [mat]() -> SparseMatrixC { return Eigen::SparseView<Eigen::MatrixXf>(mat); });
+    m.def("sparse_copy_r", [](const SparseMatrixR &m) -> SparseMatrixR { return m; });
+    m.def("sparse_copy_c", [](const SparseMatrixC &m) -> SparseMatrixC { return m; });
+    m.def("partial_copy_four_rm_r", [](const FourRowMatrixR &m) -> FourRowMatrixR { return m; });
+    m.def("partial_copy_four_rm_c", [](const FourColMatrixR &m) -> FourColMatrixR { return m; });
+    m.def("partial_copy_four_cm_r", [](const FourRowMatrixC &m) -> FourRowMatrixC { return m; });
+    m.def("partial_copy_four_cm_c", [](const FourColMatrixC &m) -> FourColMatrixC { return m; });
 
-    m.def("fixed_passthrough_r", [](const FixedMatrixR &m) -> FixedMatrixR {
+    // Test that we can cast a numpy object to a Eigen::MatrixXd explicitly
-        return m;
+    m.def("cpp_copy", [](py::handle m) { return m.cast<Eigen::MatrixXd>()(1, 0); });
-    });
+    m.def("cpp_ref_c", [](py::handle m) { return m.cast<Eigen::Ref<Eigen::MatrixXd>>()(1, 0); });
+    m.def("cpp_ref_r", [](py::handle m) { return m.cast<Eigen::Ref<MatrixXdR>>()(1, 0); });
+    m.def("cpp_ref_any", [](py::handle m) { return m.cast<py::EigenDRef<Eigen::MatrixXd>>()(1, 0); });
 
-    m.def("fixed_passthrough_c", [](const FixedMatrixC &m) -> FixedMatrixC {
-        return m;
-    });
 
-    m.def("dense_r", [mat]() -> DenseMatrixR {
+    // Test that we can prevent copying into an argument that would normally copy: First a version
-        return DenseMatrixR(mat);
+    // that would allow copying (if types or strides don't match) for comparison:
-    });
+    m.def("get_elem", &get_elem);
-
+    // Now this alternative that calls the tells pybind to fail rather than copy:
-    m.def("dense_c", [mat]() -> DenseMatrixC {
+    m.def("get_elem_nocopy", [](Eigen::Ref<const Eigen::MatrixXd> m) -> double { return get_elem(m); },
-        return DenseMatrixC(mat);
+            py::arg().noconvert());
-    });
+    // Also test a row-major-only no-copy const ref:
-
+    m.def("get_elem_rm_nocopy", [](Eigen::Ref<const Eigen::Matrix<long, -1, -1, Eigen::RowMajor>> &m) -> long { return m(2, 1); },
-    m.def("dense_passthrough_r", [](const DenseMatrixR &m) -> DenseMatrixR {
+            py::arg().noconvert());
-        return m;
-    });
-
-    m.def("dense_passthrough_c", [](const DenseMatrixC &m) -> DenseMatrixC {
-        return m;
-    });
-
-    m.def("sparse_r", [mat]() -> SparseMatrixR {
-        return Eigen::SparseView<Eigen::MatrixXf>(mat);
-    });
-
-    m.def("sparse_c", [mat]() -> SparseMatrixC {
-        return Eigen::SparseView<Eigen::MatrixXf>(mat);
-    });
-
-    m.def("sparse_passthrough_r", [](const SparseMatrixR &m) -> SparseMatrixR {
-        return m;
-    });
-
-    m.def("sparse_passthrough_c", [](const SparseMatrixC &m) -> SparseMatrixC {
-        return m;
-    });
-
-    m.def("partial_passthrough_four_rm_r", [](const FourRowMatrixR &m) -> FourRowMatrixR {
-        return m;
-    });
-    m.def("partial_passthrough_four_rm_c", [](const FourColMatrixR &m) -> FourColMatrixR {
-        return m;
-    });
-    m.def("partial_passthrough_four_cm_r", [](const FourRowMatrixC &m) -> FourRowMatrixC {
-        return m;
-    });
-    m.def("partial_passthrough_four_cm_c", [](const FourColMatrixC &m) -> FourColMatrixC {
-        return m;
-    });
 });

tests/test_eigen.py

@@ -3,7 +3,7 @@ import pytest
 with pytest.suppress(ImportError):
     import numpy as np
 
-    ref = np.array([[ 0,  3,  0,  0,  0, 11],
+    ref = np.array([[ 0.,  3,  0,  0,  0, 11],
                     [22,  0,  0,  0, 17, 11],
                     [ 7,  5,  0,  1,  0, 11],
                     [ 0,  0,  0,  0,  0, 11],
@@ -20,57 +20,130 @@ def assert_sparse_equal_ref(sparse_mat):
 
 @pytest.requires_eigen_and_numpy
 def test_fixed():
-    from pybind11_tests import fixed_r, fixed_c, fixed_passthrough_r, fixed_passthrough_c
+    from pybind11_tests import fixed_r, fixed_c, fixed_copy_r, fixed_copy_c
 
     assert_equal_ref(fixed_c())
     assert_equal_ref(fixed_r())
-    assert_equal_ref(fixed_passthrough_r(fixed_r()))
+    assert_equal_ref(fixed_copy_r(fixed_r()))
-    assert_equal_ref(fixed_passthrough_c(fixed_c()))
+    assert_equal_ref(fixed_copy_c(fixed_c()))
-    assert_equal_ref(fixed_passthrough_r(fixed_c()))
+    assert_equal_ref(fixed_copy_r(fixed_c()))
-    assert_equal_ref(fixed_passthrough_c(fixed_r()))
+    assert_equal_ref(fixed_copy_c(fixed_r()))
 
 
 @pytest.requires_eigen_and_numpy
 def test_dense():
-    from pybind11_tests import dense_r, dense_c, dense_passthrough_r, dense_passthrough_c
+    from pybind11_tests import dense_r, dense_c, dense_copy_r, dense_copy_c
 
     assert_equal_ref(dense_r())
     assert_equal_ref(dense_c())
-    assert_equal_ref(dense_passthrough_r(dense_r()))
+    assert_equal_ref(dense_copy_r(dense_r()))
-    assert_equal_ref(dense_passthrough_c(dense_c()))
+    assert_equal_ref(dense_copy_c(dense_c()))
-    assert_equal_ref(dense_passthrough_r(dense_c()))
+    assert_equal_ref(dense_copy_r(dense_c()))
-    assert_equal_ref(dense_passthrough_c(dense_r()))
+    assert_equal_ref(dense_copy_c(dense_r()))
+
 
 @pytest.requires_eigen_and_numpy
 def test_partially_fixed():
-    from pybind11_tests import partial_passthrough_four_rm_r, partial_passthrough_four_rm_c, partial_passthrough_four_cm_r, partial_passthrough_four_cm_c
+    from pybind11_tests import (partial_copy_four_rm_r, partial_copy_four_rm_c,
+                                partial_copy_four_cm_r, partial_copy_four_cm_c)
 
-    ref2 = np.array([[0,1,2,3], [4,5,6,7], [8,9,10,11], [12,13,14,15]])
+    ref2 = np.array([[0., 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11], [12, 13, 14, 15]])
-    np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2), ref2)
+    np.testing.assert_array_equal(partial_copy_four_rm_r(ref2), ref2)
-    np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2), ref2)
+    np.testing.assert_array_equal(partial_copy_four_rm_c(ref2), ref2)
-    np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2[:, 1]), ref2[:, [1]])
+    np.testing.assert_array_equal(partial_copy_four_rm_r(ref2[:, 1]), ref2[:, [1]])
-    np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2[0, :]), ref2[[0], :])
+    np.testing.assert_array_equal(partial_copy_four_rm_c(ref2[0, :]), ref2[[0], :])
-    np.testing.assert_array_equal(partial_passthrough_four_rm_r(ref2[:, (0, 2)]), ref2[:, (0,2)])
+    np.testing.assert_array_equal(partial_copy_four_rm_r(ref2[:, (0, 2)]), ref2[:, (0, 2)])
-    np.testing.assert_array_equal(partial_passthrough_four_rm_c(ref2[(3,1,2), :]), ref2[(3,1,2), :])
+    np.testing.assert_array_equal(
+        partial_copy_four_rm_c(ref2[(3, 1, 2), :]), ref2[(3, 1, 2), :])
+
+    np.testing.assert_array_equal(partial_copy_four_cm_r(ref2), ref2)
+    np.testing.assert_array_equal(partial_copy_four_cm_c(ref2), ref2)
+    np.testing.assert_array_equal(partial_copy_four_cm_r(ref2[:, 1]), ref2[:, [1]])
+    np.testing.assert_array_equal(partial_copy_four_cm_c(ref2[0, :]), ref2[[0], :])
+    np.testing.assert_array_equal(partial_copy_four_cm_r(ref2[:, (0, 2)]), ref2[:, (0, 2)])
+    np.testing.assert_array_equal(
+        partial_copy_four_cm_c(ref2[(3, 1, 2), :]), ref2[(3, 1, 2), :])
+
+
+@pytest.requires_eigen_and_numpy
+def test_mutator_descriptors():
+    from pybind11_tests import fixed_mutator_r, fixed_mutator_c, fixed_mutator_a
+    zr = np.arange(30, dtype='float32').reshape(5, 6)  # row-major
+    zc = zr.reshape(6, 5).transpose()  # column-major
+
+    fixed_mutator_r(zr)
+    fixed_mutator_c(zc)
+    fixed_mutator_a(zr)
+    fixed_mutator_a(zc)
+    with pytest.raises(TypeError) as excinfo:
+        fixed_mutator_r(zc)
+    assert ('(numpy.ndarray[float32[5, 6], flags.writeable, flags.c_contiguous]) -> arg0: None'
+            in str(excinfo.value))
+    with pytest.raises(TypeError) as excinfo:
+        fixed_mutator_c(zr)
+    assert ('(numpy.ndarray[float32[5, 6], flags.writeable, flags.f_contiguous]) -> arg0: None'
+            in str(excinfo.value))
+    with pytest.raises(TypeError) as excinfo:
+        fixed_mutator_a(np.array([[1, 2], [3, 4]], dtype='float32'))
+    assert ('(numpy.ndarray[float32[5, 6], flags.writeable]) -> arg0: None'
+            in str(excinfo.value))
+    zr.flags.writeable = False
+    with pytest.raises(TypeError):
+        fixed_mutator_r(zr)
+    with pytest.raises(TypeError):
+        fixed_mutator_a(zr)
+
+
+@pytest.requires_eigen_and_numpy
+def test_cpp_casting():
+    from pybind11_tests import (cpp_copy, cpp_ref_c, cpp_ref_r, cpp_ref_any,
+                                fixed_r, fixed_c, get_cm_ref, get_rm_ref, ReturnTester)
+    assert cpp_copy(fixed_r()) == 22.
+    assert cpp_copy(fixed_c()) == 22.
+    z = np.array([[5., 6], [7, 8]])
+    assert cpp_copy(z) == 7.
+    assert cpp_copy(get_cm_ref()) == 21.
+    assert cpp_copy(get_rm_ref()) == 21.
+    assert cpp_ref_c(get_cm_ref()) == 21.
+    assert cpp_ref_r(get_rm_ref()) == 21.
+    with pytest.raises(RuntimeError) as excinfo:
+        # Can't reference fixed_c: it contains floats, cpp_ref_any wants doubles
+        cpp_ref_any(fixed_c())
+    assert 'Unable to cast Python instance' in str(excinfo.value)
+    with pytest.raises(RuntimeError) as excinfo:
+        # Can't reference fixed_r: it contains floats, cpp_ref_any wants doubles
+        cpp_ref_any(fixed_r())
+    assert 'Unable to cast Python instance' in str(excinfo.value)
+    assert cpp_ref_any(ReturnTester.create()) == 1.
+
+    assert cpp_ref_any(get_cm_ref()) == 21.
+    assert cpp_ref_any(get_cm_ref()) == 21.
+
+
+@pytest.requires_eigen_and_numpy
+def test_pass_readonly_array():
+    from pybind11_tests import fixed_copy_r, fixed_r, fixed_r_const
+    z = np.full((5, 6), 42.0)
+    z.flags.writeable = False
+    np.testing.assert_array_equal(z, fixed_copy_r(z))
+    np.testing.assert_array_equal(fixed_r_const(), fixed_r())
+    assert not fixed_r_const().flags.writeable
+    np.testing.assert_array_equal(fixed_copy_r(fixed_r_const()), fixed_r_const())
 
-    np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2), ref2)
-    np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2), ref2)
-    np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2[:, 1]), ref2[:, [1]])
-    np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2[0, :]), ref2[[0], :])
-    np.testing.assert_array_equal(partial_passthrough_four_cm_r(ref2[:, (0, 2)]), ref2[:, (0,2)])
-    np.testing.assert_array_equal(partial_passthrough_four_cm_c(ref2[(3,1,2), :]), ref2[(3,1,2), :])
 
 @pytest.requires_eigen_and_numpy
 def test_nonunit_stride_from_python():
-    from pybind11_tests import double_row, double_col, double_mat_cm, double_mat_rm
+    from pybind11_tests import (
+        double_row, double_col, double_mat_cm, double_mat_rm,
+        double_threec, double_threer)
 
     counting_mat = np.arange(9.0, dtype=np.float32).reshape((3, 3))
-    first_row = counting_mat[0, :]
+    second_row = counting_mat[1, :]
-    first_col = counting_mat[:, 0]
+    second_col = counting_mat[:, 1]
-    np.testing.assert_array_equal(double_row(first_row), 2.0 * first_row)
+    np.testing.assert_array_equal(double_row(second_row), 2.0 * second_row)
-    np.testing.assert_array_equal(double_col(first_row), 2.0 * first_row)
+    np.testing.assert_array_equal(double_col(second_row), 2.0 * second_row)
-    np.testing.assert_array_equal(double_row(first_col), 2.0 * first_col)
+    np.testing.assert_array_equal(double_row(second_col), 2.0 * second_col)
-    np.testing.assert_array_equal(double_col(first_col), 2.0 * first_col)
+    np.testing.assert_array_equal(double_col(second_col), 2.0 * second_col)
 
     counting_3d = np.arange(27.0, dtype=np.float32).reshape((3, 3, 3))
     slices = [counting_3d[0, :, :], counting_3d[:, 0, :], counting_3d[:, :, 0]]
@@ -78,6 +151,11 @@ def test_nonunit_stride_from_python():
         np.testing.assert_array_equal(double_mat_cm(ref_mat), 2.0 * ref_mat)
         np.testing.assert_array_equal(double_mat_rm(ref_mat), 2.0 * ref_mat)
 
+    # Mutator:
+    double_threer(second_row)
+    double_threec(second_col)
+    np.testing.assert_array_equal(counting_mat, [[0., 2, 2], [6, 16, 10], [6, 14, 8]])
+
 
 @pytest.requires_eigen_and_numpy
 def test_nonunit_stride_to_python():
@@ -95,21 +173,396 @@ def test_nonunit_stride_to_python():
 
 @pytest.requires_eigen_and_numpy
 def test_eigen_ref_to_python():
-    from pybind11_tests import cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6
+    from pybind11_tests import cholesky1, cholesky2, cholesky3, cholesky4
 
-    chols = [cholesky1, cholesky2, cholesky3, cholesky4, cholesky5, cholesky6]
+    chols = [cholesky1, cholesky2, cholesky3, cholesky4]
     for i, chol in enumerate(chols, start=1):
-        mymat = chol(np.array([[1, 2, 4], [2, 13, 23], [4, 23, 77]]))
+        mymat = chol(np.array([[1., 2, 4], [2, 13, 23], [4, 23, 77]]))
         assert np.all(mymat == np.array([[1, 0, 0], [2, 3, 0], [4, 5, 6]])), "cholesky{}".format(i)
 
 
+def assign_both(a1, a2, r, c, v):
+    a1[r, c] = v
+    a2[r, c] = v
+
+
+def array_copy_but_one(a, r, c, v):
+    z = np.array(a, copy=True)
+    z[r, c] = v
+    return z
+
+
+@pytest.requires_eigen_and_numpy
+def test_eigen_return_references():
+    """Tests various ways of returning references and non-referencing copies"""
+    from pybind11_tests import ReturnTester
+    master = np.ones((10, 10))
+    a = ReturnTester()
+    a_get1 = a.get()
+    assert not a_get1.flags.owndata and a_get1.flags.writeable
+    assign_both(a_get1, master, 3, 3, 5)
+    a_get2 = a.get_ptr()
+    assert not a_get2.flags.owndata and a_get2.flags.writeable
+    assign_both(a_get1, master, 2, 3, 6)
+
+    a_view1 = a.view()
+    assert not a_view1.flags.owndata and not a_view1.flags.writeable
+    with pytest.raises(ValueError):
+        a_view1[2, 3] = 4
+    a_view2 = a.view_ptr()
+    assert not a_view2.flags.owndata and not a_view2.flags.writeable
+    with pytest.raises(ValueError):
+        a_view2[2, 3] = 4
+
+    a_copy1 = a.copy_get()
+    assert a_copy1.flags.owndata and a_copy1.flags.writeable
+    np.testing.assert_array_equal(a_copy1, master)
+    a_copy1[7, 7] = -44  # Shouldn't affect anything else
+    c1want = array_copy_but_one(master, 7, 7, -44)
+    a_copy2 = a.copy_view()
+    assert a_copy2.flags.owndata and a_copy2.flags.writeable
+    np.testing.assert_array_equal(a_copy2, master)
+    a_copy2[4, 4] = -22  # Shouldn't affect anything else
+    c2want = array_copy_but_one(master, 4, 4, -22)
+
+    a_ref1 = a.ref()
+    assert not a_ref1.flags.owndata and a_ref1.flags.writeable
+    assign_both(a_ref1, master, 1, 1, 15)
+    a_ref2 = a.ref_const()
+    assert not a_ref2.flags.owndata and not a_ref2.flags.writeable
+    with pytest.raises(ValueError):
+        a_ref2[5, 5] = 33
+    a_ref3 = a.ref_safe()
+    assert not a_ref3.flags.owndata and a_ref3.flags.writeable
+    assign_both(a_ref3, master, 0, 7, 99)
+    a_ref4 = a.ref_const_safe()
+    assert not a_ref4.flags.owndata and not a_ref4.flags.writeable
+    with pytest.raises(ValueError):
+        a_ref4[7, 0] = 987654321
+
+    a_copy3 = a.copy_ref()
+    assert a_copy3.flags.owndata and a_copy3.flags.writeable
+    np.testing.assert_array_equal(a_copy3, master)
+    a_copy3[8, 1] = 11
+    c3want = array_copy_but_one(master, 8, 1, 11)
+    a_copy4 = a.copy_ref_const()
+    assert a_copy4.flags.owndata and a_copy4.flags.writeable
+    np.testing.assert_array_equal(a_copy4, master)
+    a_copy4[8, 4] = 88
+    c4want = array_copy_but_one(master, 8, 4, 88)
+
+    a_block1 = a.block(3, 3, 2, 2)
+    assert not a_block1.flags.owndata and a_block1.flags.writeable
+    a_block1[0, 0] = 55
+    master[3, 3] = 55
+    a_block2 = a.block_safe(2, 2, 3, 2)
+    assert not a_block2.flags.owndata and a_block2.flags.writeable
+    a_block2[2, 1] = -123
+    master[4, 3] = -123
+    a_block3 = a.block_const(6, 7, 4, 3)
+    assert not a_block3.flags.owndata and not a_block3.flags.writeable
+    with pytest.raises(ValueError):
+        a_block3[2, 2] = -44444
+
+    a_copy5 = a.copy_block(2, 2, 2, 3)
+    assert a_copy5.flags.owndata and a_copy5.flags.writeable
+    np.testing.assert_array_equal(a_copy5, master[2:4, 2:5])
+    a_copy5[1, 1] = 777
+    c5want = array_copy_but_one(master[2:4, 2:5], 1, 1, 777)
+
+    a_corn1 = a.corners()
+    assert not a_corn1.flags.owndata and a_corn1.flags.writeable
+    a_corn1 *= 50
+    a_corn1[1, 1] = 999
+    master[0, 0] = 50
+    master[0, 9] = 50
+    master[9, 0] = 50
+    master[9, 9] = 999
+    a_corn2 = a.corners_const()
+    assert not a_corn2.flags.owndata and not a_corn2.flags.writeable
+    with pytest.raises(ValueError):
+        a_corn2[1, 0] = 51
+
+    # All of the changes made all the way along should be visible everywhere
+    # now (except for the copies, of course)
+    np.testing.assert_array_equal(a_get1, master)
+    np.testing.assert_array_equal(a_get2, master)
+    np.testing.assert_array_equal(a_view1, master)
+    np.testing.assert_array_equal(a_view2, master)
+    np.testing.assert_array_equal(a_ref1, master)
+    np.testing.assert_array_equal(a_ref2, master)
+    np.testing.assert_array_equal(a_ref3, master)
+    np.testing.assert_array_equal(a_ref4, master)
+    np.testing.assert_array_equal(a_block1, master[3:5, 3:5])
+    np.testing.assert_array_equal(a_block2, master[2:5, 2:4])
+    np.testing.assert_array_equal(a_block3, master[6:10, 7:10])
+    np.testing.assert_array_equal(a_corn1, master[0::master.shape[0] - 1, 0::master.shape[1] - 1])
+    np.testing.assert_array_equal(a_corn2, master[0::master.shape[0] - 1, 0::master.shape[1] - 1])
+
+    np.testing.assert_array_equal(a_copy1, c1want)
+    np.testing.assert_array_equal(a_copy2, c2want)
+    np.testing.assert_array_equal(a_copy3, c3want)
+    np.testing.assert_array_equal(a_copy4, c4want)
+    np.testing.assert_array_equal(a_copy5, c5want)
+
+
+def assert_keeps_alive(cl, method, *args):
+    from pybind11_tests import ConstructorStats
+    cstats = ConstructorStats.get(cl)
+    start_with = cstats.alive()
+    a = cl()
+    assert cstats.alive() == start_with + 1
+    z = method(a, *args)
+    assert cstats.alive() == start_with + 1
+    del a
+    # Here's the keep alive in action:
+    assert cstats.alive() == start_with + 1
+    del z
+    # Keep alive should have expired:
+    assert cstats.alive() == start_with
+
+
+@pytest.requires_eigen_and_numpy
+def test_eigen_keepalive():
+    from pybind11_tests import ReturnTester, ConstructorStats
+    a = ReturnTester()
+
+    cstats = ConstructorStats.get(ReturnTester)
+    assert cstats.alive() == 1
+    unsafe = [a.ref(), a.ref_const(), a.block(1, 2, 3, 4)]
+    copies = [a.copy_get(), a.copy_view(), a.copy_ref(), a.copy_ref_const(),
+              a.copy_block(4, 3, 2, 1)]
+    del a
+    assert cstats.alive() == 0
+    del unsafe
+    del copies
+
+    for meth in [ReturnTester.get, ReturnTester.get_ptr, ReturnTester.view,
+                 ReturnTester.view_ptr, ReturnTester.ref_safe, ReturnTester.ref_const_safe,
+                 ReturnTester.corners, ReturnTester.corners_const]:
+        assert_keeps_alive(ReturnTester, meth)
+
+    for meth in [ReturnTester.block_safe, ReturnTester.block_const]:
+        assert_keeps_alive(ReturnTester, meth, 4, 3, 2, 1)
+
+
+@pytest.requires_eigen_and_numpy
+def test_eigen_ref_mutators():
+    """Tests whether Eigen can mutate numpy values"""
+    from pybind11_tests import add_rm, add_cm, add_any, add1, add2
+    orig = np.array([[1., 2, 3], [4, 5, 6], [7, 8, 9]])
+    zr = np.array(orig)
+    zc = np.array(orig, order='F')
+    add_rm(zr, 1, 0, 100)
+    assert np.all(zr == np.array([[1., 2, 3], [104, 5, 6], [7, 8, 9]]))
+    add_cm(zc, 1, 0, 200)
+    assert np.all(zc == np.array([[1., 2, 3], [204, 5, 6], [7, 8, 9]]))
+
+    add_any(zr, 1, 0, 20)
+    assert np.all(zr == np.array([[1., 2, 3], [124, 5, 6], [7, 8, 9]]))
+    add_any(zc, 1, 0, 10)
+    assert np.all(zc == np.array([[1., 2, 3], [214, 5, 6], [7, 8, 9]]))
+
+    # Can't reference a col-major array with a row-major Ref, and vice versa:
+    with pytest.raises(TypeError):
+        add_rm(zc, 1, 0, 1)
+    with pytest.raises(TypeError):
+        add_cm(zr, 1, 0, 1)
+
+    # Overloads:
+    add1(zr, 1, 0, -100)
+    add2(zr, 1, 0, -20)
+    assert np.all(zr == orig)
+    add1(zc, 1, 0, -200)
+    add2(zc, 1, 0, -10)
+    assert np.all(zc == orig)
+
+    # a non-contiguous slice (this won't work on either the row- or
+    # column-contiguous refs, but should work for the any)
+    cornersr = zr[0::2, 0::2]
+    cornersc = zc[0::2, 0::2]
+
+    assert np.all(cornersr == np.array([[1., 3], [7, 9]]))
+    assert np.all(cornersc == np.array([[1., 3], [7, 9]]))
+
+    with pytest.raises(TypeError):
+        add_rm(cornersr, 0, 1, 25)
+    with pytest.raises(TypeError):
+        add_cm(cornersr, 0, 1, 25)
+    with pytest.raises(TypeError):
+        add_rm(cornersc, 0, 1, 25)
+    with pytest.raises(TypeError):
+        add_cm(cornersc, 0, 1, 25)
+    add_any(cornersr, 0, 1, 25)
+    add_any(cornersc, 0, 1, 44)
+    assert np.all(zr == np.array([[1., 2, 28], [4, 5, 6], [7, 8, 9]]))
+    assert np.all(zc == np.array([[1., 2, 47], [4, 5, 6], [7, 8, 9]]))
+
+    # You shouldn't be allowed to pass a non-writeable array to a mutating Eigen method:
+    zro = zr[0:4, 0:4]
+    zro.flags.writeable = False
+    with pytest.raises(TypeError):
+        add_rm(zro, 0, 0, 0)
+    with pytest.raises(TypeError):
+        add_any(zro, 0, 0, 0)
+    with pytest.raises(TypeError):
+        add1(zro, 0, 0, 0)
+    with pytest.raises(TypeError):
+        add2(zro, 0, 0, 0)
+
+    # integer array shouldn't be passable to a double-matrix-accepting mutating func:
+    zi = np.array([[1, 2], [3, 4]])
+    with pytest.raises(TypeError):
+        add_rm(zi)
+
+
+@pytest.requires_eigen_and_numpy
+def test_numpy_ref_mutators():
+    """Tests numpy mutating Eigen matrices (for returned Eigen::Ref<...>s)"""
+    from pybind11_tests import (
+        get_cm_ref, get_cm_const_ref, get_rm_ref, get_rm_const_ref, reset_refs)
+    reset_refs()  # In case another test already changed it
+
+    zc = get_cm_ref()
+    zcro = get_cm_const_ref()
+    zr = get_rm_ref()
+    zrro = get_rm_const_ref()
+
+    assert [zc[1, 2], zcro[1, 2], zr[1, 2], zrro[1, 2]] == [23] * 4
+
+    assert not zc.flags.owndata and zc.flags.writeable
+    assert not zr.flags.owndata and zr.flags.writeable
+    assert not zcro.flags.owndata and not zcro.flags.writeable
+    assert not zrro.flags.owndata and not zrro.flags.writeable
+
+    zc[1, 2] = 99
+    expect = np.array([[11., 12, 13], [21, 22, 99], [31, 32, 33]])
+    # We should have just changed zc, of course, but also zcro and the original eigen matrix
+    assert np.all(zc == expect)
+    assert np.all(zcro == expect)
+    assert np.all(get_cm_ref() == expect)
+
+    zr[1, 2] = 99
+    assert np.all(zr == expect)
+    assert np.all(zrro == expect)
+    assert np.all(get_rm_ref() == expect)
+
+    # Make sure the readonly ones are numpy-readonly:
+    with pytest.raises(ValueError):
+        zcro[1, 2] = 6
+    with pytest.raises(ValueError):
+        zrro[1, 2] = 6
+
+    # We should be able to explicitly copy like this (and since we're copying,
+    # the const should drop away)
+    y1 = np.array(get_cm_const_ref())
+
+    assert y1.flags.owndata and y1.flags.writeable
+    # We should get copies of the eigen data, which was modified above:
+    assert y1[1, 2] == 99
+    y1[1, 2] += 12
+    assert y1[1, 2] == 111
+    assert zc[1, 2] == 99  # Make sure we aren't referencing the original
+
+
+@pytest.requires_eigen_and_numpy
+def test_both_ref_mutators():
+    """Tests a complex chain of nested eigen/numpy references"""
+    from pybind11_tests import (
+        incr_matrix, get_cm_ref, incr_matrix_any, even_cols, even_rows, reset_refs)
+    reset_refs()  # In case another test already changed it
+
+    z = get_cm_ref()  # numpy -> eigen
+    z[0, 2] -= 3
+    z2 = incr_matrix(z, 1)  # numpy -> eigen -> numpy -> eigen
+    z2[1, 1] += 6
+    z3 = incr_matrix(z, 2)  # (numpy -> eigen)^3
+    z3[2, 2] += -5
+    z4 = incr_matrix(z, 3)  # (numpy -> eigen)^4
+    z4[1, 1] -= 1
+    z5 = incr_matrix(z, 4)  # (numpy -> eigen)^5
+    z5[0, 0] = 0
+    assert np.all(z == z2)
+    assert np.all(z == z3)
+    assert np.all(z == z4)
+    assert np.all(z == z5)
+    expect = np.array([[0., 22, 20], [31, 37, 33], [41, 42, 38]])
+    assert np.all(z == expect)
+
+    y = np.array(range(100), dtype='float64').reshape(10, 10)
+    y2 = incr_matrix_any(y, 10)  # np -> eigen -> np
+    y3 = incr_matrix_any(y2[0::2, 0::2], -33)  # np -> eigen -> np slice -> np -> eigen -> np
+    y4 = even_rows(y3)  # numpy -> eigen slice -> (... y3)
+    y5 = even_cols(y4)  # numpy -> eigen slice -> (... y4)
+    y6 = incr_matrix_any(y5, 1000)  # numpy -> eigen -> (... y5)
+
+    # Apply same mutations using just numpy:
+    yexpect = np.array(range(100), dtype='float64').reshape(10, 10)
+    yexpect += 10
+    yexpect[0::2, 0::2] -= 33
+    yexpect[0::4, 0::4] += 1000
+    assert np.all(y6 == yexpect[0::4, 0::4])
+    assert np.all(y5 == yexpect[0::4, 0::4])
+    assert np.all(y4 == yexpect[0::4, 0::2])
+    assert np.all(y3 == yexpect[0::2, 0::2])
+    assert np.all(y2 == yexpect)
+    assert np.all(y == yexpect)
+
+
+@pytest.requires_eigen_and_numpy
+def test_nocopy_wrapper():
+    from pybind11_tests import get_elem, get_elem_nocopy, get_elem_rm_nocopy
+    # get_elem requires a column-contiguous matrix reference, but should be
+    # callable with other types of matrix (via copying):
+    int_matrix_colmajor = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], order='F')
+    dbl_matrix_colmajor = np.array(int_matrix_colmajor, dtype='double', order='F', copy=True)
+    int_matrix_rowmajor = np.array(int_matrix_colmajor, order='C', copy=True)
+    dbl_matrix_rowmajor = np.array(int_matrix_rowmajor, dtype='double', order='C', copy=True)
+
+    # All should be callable via get_elem:
+    assert get_elem(int_matrix_colmajor) == 8
+    assert get_elem(dbl_matrix_colmajor) == 8
+    assert get_elem(int_matrix_rowmajor) == 8
+    assert get_elem(dbl_matrix_rowmajor) == 8
+
+    # All but the second should fail with get_elem_nocopy:
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_nocopy(int_matrix_colmajor)
+    assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.f_contiguous' in str(excinfo.value))
+    assert get_elem_nocopy(dbl_matrix_colmajor) == 8
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_nocopy(int_matrix_rowmajor)
+    assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.f_contiguous' in str(excinfo.value))
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_nocopy(dbl_matrix_rowmajor)
+    assert ('get_elem_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.f_contiguous' in str(excinfo.value))
+
+    # For the row-major test, we take a long matrix in row-major, so only the third is allowed:
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_rm_nocopy(int_matrix_colmajor)
+    assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.c_contiguous' in str(excinfo.value))
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_rm_nocopy(dbl_matrix_colmajor)
+    assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.c_contiguous' in str(excinfo.value))
+    assert get_elem_rm_nocopy(int_matrix_rowmajor) == 8
+    with pytest.raises(TypeError) as excinfo:
+        get_elem_rm_nocopy(dbl_matrix_rowmajor)
+    assert ('get_elem_rm_nocopy(): incompatible function arguments.' in str(excinfo.value) and
+            ', flags.c_contiguous' in str(excinfo.value))
+
+
 @pytest.requires_eigen_and_numpy
 def test_special_matrix_objects():
     from pybind11_tests import incr_diag, symmetric_upper, symmetric_lower
 
-    assert np.all(incr_diag(7) == np.diag([1, 2, 3, 4, 5, 6, 7]))
+    assert np.all(incr_diag(7) == np.diag([1., 2, 3, 4, 5, 6, 7]))
 
-    asymm = np.array([[ 1,  2,  3,  4],
+    asymm = np.array([[ 1.,  2,  3,  4],
                       [ 5,  6,  7,  8],
                       [ 9, 10, 11, 12],
                       [13, 14, 15, 16]])
@@ -141,23 +594,23 @@ def test_dense_signature(doc):
 
 @pytest.requires_eigen_and_scipy
 def test_sparse():
-    from pybind11_tests import sparse_r, sparse_c, sparse_passthrough_r, sparse_passthrough_c
+    from pybind11_tests import sparse_r, sparse_c, sparse_copy_r, sparse_copy_c
 
     assert_sparse_equal_ref(sparse_r())
     assert_sparse_equal_ref(sparse_c())
-    assert_sparse_equal_ref(sparse_passthrough_r(sparse_r()))
+    assert_sparse_equal_ref(sparse_copy_r(sparse_r()))
-    assert_sparse_equal_ref(sparse_passthrough_c(sparse_c()))
+    assert_sparse_equal_ref(sparse_copy_c(sparse_c()))
-    assert_sparse_equal_ref(sparse_passthrough_r(sparse_c()))
+    assert_sparse_equal_ref(sparse_copy_r(sparse_c()))
-    assert_sparse_equal_ref(sparse_passthrough_c(sparse_r()))
+    assert_sparse_equal_ref(sparse_copy_c(sparse_r()))
 
 
 @pytest.requires_eigen_and_scipy
 def test_sparse_signature(doc):
-    from pybind11_tests import sparse_passthrough_r, sparse_passthrough_c
+    from pybind11_tests import sparse_copy_r, sparse_copy_c
 
-    assert doc(sparse_passthrough_r) == """
+    assert doc(sparse_copy_r) == """
-        sparse_passthrough_r(arg0: scipy.sparse.csr_matrix[float32]) -> scipy.sparse.csr_matrix[float32]
+        sparse_copy_r(arg0: scipy.sparse.csr_matrix[float32]) -> scipy.sparse.csr_matrix[float32]
     """  # noqa: E501 line too long
-    assert doc(sparse_passthrough_c) == """
+    assert doc(sparse_copy_c) == """
-        sparse_passthrough_c(arg0: scipy.sparse.csc_matrix[float32]) -> scipy.sparse.csc_matrix[float32]
+        sparse_copy_c(arg0: scipy.sparse.csc_matrix[float32]) -> scipy.sparse.csc_matrix[float32]
     """  # noqa: E501 line too long