unittest/bind_virtual_factory.cpp

+/// Copyright 2023 LAAS-CNRS, INRIA
+#include <eigenpy/eigenpy.hpp>
+
+using std::shared_ptr;
+namespace bp = boost::python;
+
+// fwd declaration
+struct MyVirtualData;
+
+/// A virtual class with two pure virtual functions taking different signatures,
+/// and a polymorphic factory function.
+struct MyVirtualClass {
+  MyVirtualClass() {}
+  virtual ~MyVirtualClass() {}
+
+  // polymorphic fn taking arg by shared_ptr
+  virtual void doSomethingPtr(shared_ptr<MyVirtualData> const &data) const = 0;
+  // polymorphic fn taking arg by reference
+  virtual void doSomethingRef(MyVirtualData &data) const = 0;
+
+  virtual shared_ptr<MyVirtualData> createData() const {
+    return std::make_shared<MyVirtualData>(*this);
+  }
+};
+
+struct MyVirtualData {
+  MyVirtualData(MyVirtualClass const &) {}
+  virtual ~MyVirtualData() {}  // virtual dtor to mark class as polymorphic
+};
+
+shared_ptr<MyVirtualData> callDoSomethingPtr(const MyVirtualClass &obj) {
+  auto d = obj.createData();
+  printf("Created MyVirtualData with address %p\n", (void *)d.get());
+  obj.doSomethingPtr(d);
+  return d;
+}
+
+shared_ptr<MyVirtualData> callDoSomethingRef(const MyVirtualClass &obj) {
+  auto d = obj.createData();
+  printf("Created MyVirtualData with address %p\n", (void *)d.get());
+  obj.doSomethingRef(*d);
+  return d;
+}
+
+void throw_virtual_not_implemented_error() {
+  throw std::runtime_error("Called C++ virtual function.");
+}
+
+/// Wrapper classes
+struct VirtualClassWrapper : MyVirtualClass, bp::wrapper<MyVirtualClass> {
+  void doSomethingPtr(shared_ptr<MyVirtualData> const &data) const override {
+    if (bp::override fo = this->get_override("doSomethingPtr")) {
+      /// shared_ptr HAS to be passed by value.
+      /// Boost.Python's argument converter has the wrong behaviour for
+      /// reference_wrapper<shared_ptr<T>>, so boost::ref(data) does not work.
+      fo(data);
+      return;
+    }
+    throw_virtual_not_implemented_error();
+  }
+
+  /// The data object is passed by mutable reference to this function,
+  /// and wrapped in a @c boost::reference_wrapper when passed to the override.
+  /// Otherwise, Boost.Python's argument converter will convert to Python by
+  /// value and create a copy.
+  void doSomethingRef(MyVirtualData &data) const override {
+    if (bp::override fo = this->get_override("doSomethingRef")) {
+      fo(boost::ref(data));
+      return;
+    }
+    throw_virtual_not_implemented_error();
+  }
+
+  shared_ptr<MyVirtualData> createData() const override {
+    if (bp::override fo = this->get_override("createData")) {
+      bp::object result = fo().as<bp::object>();
+      return bp::extract<shared_ptr<MyVirtualData> >(result);
+    }
+    return default_createData();
+  }
+
+  shared_ptr<MyVirtualData> default_createData() const {
+    return MyVirtualClass::createData();
+  }
+};
+
+/// This "trampoline class" does nothing but is ABSOLUTELY required to ensure
+/// downcasting works properly with non-smart ptr signatures. Otherwise,
+/// there is no handle to the original Python object ( @c PyObject *).
+/// Every single polymorphic type exposed to Python should be exposed through
+/// such a trampoline. Users can also create their own wrapper classes by taking
+/// inspiration from boost::python::wrapper<T>.
+struct DataWrapper : MyVirtualData, bp::wrapper<MyVirtualData> {
+  /// we have to use-declare non-defaulted constructors
+  /// (see https://en.cppreference.com/w/cpp/language/default_constructor)
+  /// or define them manually.
+  using MyVirtualData::MyVirtualData;
+};
+
+/// Take and return a const reference
+const MyVirtualData &iden_ref(const MyVirtualData &d) {
+  // try cast to holder
+  return d;
+}
+
+/// Take a shared_ptr (by const reference or value, doesn't matter), return by
+/// const reference
+const MyVirtualData &iden_shared(const shared_ptr<MyVirtualData> &d) {
+  // get boost.python's custom deleter
+  // boost.python hides the handle to the original object in there
+  // dter being nonzero indicates shared_ptr was wrapped by Boost.Python
+  auto *dter = std::get_deleter<bp::converter::shared_ptr_deleter>(d);
+  if (dter != 0) printf("> input shared_ptr has a deleter\n");
+  return *d;
+}
+
+/// Take and return a shared_ptr
+shared_ptr<MyVirtualData> copy_shared(const shared_ptr<MyVirtualData> &d) {
+  auto *dter = std::get_deleter<bp::converter::shared_ptr_deleter>(d);
+  if (dter != 0) printf("> input shared_ptr has a deleter\n");
+  return d;
+}
+
+BOOST_PYTHON_MODULE(bind_virtual_factory) {
+  assert(std::is_polymorphic<MyVirtualClass>::value &&
+         "MyVirtualClass should be polymorphic!");
+  assert(std::is_polymorphic<MyVirtualData>::value &&
+         "MyVirtualData should be polymorphic!");
+
+  bp::class_<VirtualClassWrapper, boost::noncopyable>(
+      "MyVirtualClass", bp::init<>(bp::args("self")))
+      .def("doSomething", bp::pure_virtual(&MyVirtualClass::doSomethingPtr),
+           bp::args("self", "data"))
+      .def("doSomethingRef", bp::pure_virtual(&MyVirtualClass::doSomethingRef),
+           bp::args("self", "data"))
+      .def("createData", &MyVirtualClass::createData,
+           &VirtualClassWrapper::default_createData, bp::args("self"));
+
+  bp::register_ptr_to_python<shared_ptr<MyVirtualData> >();
+  /// Trampoline used as 1st argument
+  /// otherwise if passed as "HeldType", we need to define
+  /// the constructor and call initializer manually.
+  bp::class_<DataWrapper, boost::noncopyable>("MyVirtualData", bp::no_init)
+      .def(bp::init<MyVirtualClass const &>(bp::args("self", "model")));
+
+  bp::def("callDoSomethingPtr", callDoSomethingPtr, bp::args("obj"));
+  bp::def("callDoSomethingRef", callDoSomethingRef, bp::args("obj"));
+
+  bp::def("iden_ref", iden_ref, bp::return_internal_reference<>());
+  bp::def("iden_shared", iden_shared, bp::return_internal_reference<>());
+  bp::def("copy_shared", copy_shared);
+}

unittest/complex.cpp

@@ -4,76 +4,84 @@
 
 #include "eigenpy/eigenpy.hpp"
 
-namespace Eigen
-{
-  #define EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix)   \
-  /** \ingroup matrixtypedefs */                                    \
-  typedef Matrix<Type, Size, Size> Matrix##SizeSuffix##TypeSuffix;  \
-  /** \ingroup matrixtypedefs */                                    \
-  typedef Matrix<Type, Size, 1>    Vector##SizeSuffix##TypeSuffix;  \
-  /** \ingroup matrixtypedefs */                                    \
-  typedef Matrix<Type, 1, Size>    RowVector##SizeSuffix##TypeSuffix;
+namespace Eigen {
+#define EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix)    \
+  /** \ingroup matrixtypedefs */                                   \
+  typedef Matrix<Type, Size, Size> Matrix##SizeSuffix##TypeSuffix; \
+  /** \ingroup matrixtypedefs */                                   \
+  typedef Matrix<Type, Size, 1> Vector##SizeSuffix##TypeSuffix;    \
+  /** \ingroup matrixtypedefs */                                   \
+  typedef Matrix<Type, 1, Size> RowVector##SizeSuffix##TypeSuffix;
 
-  #define EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, Size)         \
-  /** \ingroup matrixtypedefs */                                    \
-  typedef Matrix<Type, Size, Dynamic> Matrix##Size##X##TypeSuffix;  \
-  /** \ingroup matrixtypedefs */                                    \
+#define EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, Size)          \
+  /** \ingroup matrixtypedefs */                                   \
+  typedef Matrix<Type, Size, Dynamic> Matrix##Size##X##TypeSuffix; \
+  /** \ingroup matrixtypedefs */                                   \
   typedef Matrix<Type, Dynamic, Size> Matrix##X##Size##TypeSuffix;
 
-  #define EIGEN_MAKE_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
-  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 2, 2) \
-  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 3, 3) \
-  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 4, 4) \
-  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Dynamic, X) \
-  EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 2) \
-  EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 3) \
+#define EIGEN_MAKE_TYPEDEFS_ALL_SIZES(Type, TypeSuffix) \
+  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 2, 2)           \
+  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 3, 3)           \
+  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, 4, 4)           \
+  EIGEN_MAKE_TYPEDEFS(Type, TypeSuffix, Dynamic, X)     \
+  EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 2)        \
+  EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 3)        \
   EIGEN_MAKE_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
 
-  EIGEN_MAKE_TYPEDEFS_ALL_SIZES(long double,               ld)
-  EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<long double>, cld)
+EIGEN_MAKE_TYPEDEFS_ALL_SIZES(long double, ld)
+EIGEN_MAKE_TYPEDEFS_ALL_SIZES(std::complex<long double>, cld)
 
-  #undef EIGEN_MAKE_TYPEDEFS_ALL_SIZES
-  #undef EIGEN_MAKE_TYPEDEFS
-  #undef EIGEN_MAKE_FIXED_TYPEDEFS
-}
+#undef EIGEN_MAKE_TYPEDEFS_ALL_SIZES
+#undef EIGEN_MAKE_TYPEDEFS
+#undef EIGEN_MAKE_FIXED_TYPEDEFS
+}  // namespace Eigen
 
-template<typename ComplexMatrix>
-typename Eigen::Matrix<typename ComplexMatrix::RealScalar,ComplexMatrix::RowsAtCompileTime,ComplexMatrix::ColsAtCompileTime,ComplexMatrix::Options>
-real(const Eigen::MatrixBase<ComplexMatrix> & complex_mat)
-{
+template <typename ComplexMatrix>
+typename Eigen::Matrix<typename ComplexMatrix::RealScalar,
+                       ComplexMatrix::RowsAtCompileTime,
+                       ComplexMatrix::ColsAtCompileTime, ComplexMatrix::Options>
+real(const Eigen::MatrixBase<ComplexMatrix> &complex_mat) {
   return complex_mat.real();
 }
 
-template<typename ComplexMatrix>
-typename Eigen::Matrix<typename ComplexMatrix::RealScalar,ComplexMatrix::RowsAtCompileTime,ComplexMatrix::ColsAtCompileTime,ComplexMatrix::Options>
-imag(const Eigen::MatrixBase<ComplexMatrix> & complex_mat)
-{
+template <typename ComplexMatrix>
+typename Eigen::Matrix<typename ComplexMatrix::RealScalar,
+                       ComplexMatrix::RowsAtCompileTime,
+                       ComplexMatrix::ColsAtCompileTime, ComplexMatrix::Options>
+imag(const Eigen::MatrixBase<ComplexMatrix> &complex_mat) {
   return complex_mat.imag();
 }
 
-template<typename Scalar, int Rows, int Cols, int Options>
-Eigen::Matrix<std::complex<Scalar>,Rows,Cols,Options>
-ascomplex(const Eigen::Matrix<Scalar,Rows,Cols,Options> & mat)
-{
-  typedef Eigen::Matrix<std::complex<Scalar>,Rows,Cols,Options> ReturnType;
+template <typename Scalar, int Rows, int Cols, int Options>
+Eigen::Matrix<std::complex<Scalar>, Rows, Cols, Options> ascomplex(
+    const Eigen::Matrix<Scalar, Rows, Cols, Options> &mat) {
+  typedef Eigen::Matrix<std::complex<Scalar>, Rows, Cols, Options> ReturnType;
   return ReturnType(mat.template cast<std::complex<Scalar> >());
 }
 
-BOOST_PYTHON_MODULE(complex)
-{
+BOOST_PYTHON_MODULE(complex) {
   using namespace Eigen;
   namespace bp = boost::python;
   eigenpy::enableEigenPy();
-  
-  bp::def("ascomplex", ascomplex<float,Eigen::Dynamic,Eigen::Dynamic,0>);
-  bp::def("ascomplex", ascomplex<double,Eigen::Dynamic,Eigen::Dynamic,0>);
-  bp::def("ascomplex", ascomplex<long double,Eigen::Dynamic,Eigen::Dynamic,0>);
-  
-  bp::def("real", (MatrixXf (*)(const Eigen::MatrixBase<MatrixXcf> &))&real<MatrixXcf>);
-  bp::def("real", (MatrixXd (*)(const Eigen::MatrixBase<MatrixXcd> &))&real<MatrixXcd>);
-  bp::def("real", (MatrixXld (*)(const Eigen::MatrixBase<MatrixXcld> &))&real<MatrixXcld>);
-  
-  bp::def("imag", (MatrixXf (*)(const Eigen::MatrixBase<MatrixXcf> &))&imag<MatrixXcf>);
-  bp::def("imag", (MatrixXd (*)(const Eigen::MatrixBase<MatrixXcd> &))&imag<MatrixXcd>);
-  bp::def("imag", (MatrixXld (*)(const Eigen::MatrixBase<MatrixXcld> &))&imag<MatrixXcld>);
+
+  bp::def("ascomplex", ascomplex<float, Eigen::Dynamic, Eigen::Dynamic, 0>);
+  bp::def("ascomplex", ascomplex<double, Eigen::Dynamic, Eigen::Dynamic, 0>);
+  bp::def("ascomplex",
+          ascomplex<long double, Eigen::Dynamic, Eigen::Dynamic, 0>);
+
+  bp::def("real",
+          (MatrixXf(*)(const Eigen::MatrixBase<MatrixXcf> &))&real<MatrixXcf>);
+  bp::def("real",
+          (MatrixXd(*)(const Eigen::MatrixBase<MatrixXcd> &))&real<MatrixXcd>);
+  bp::def(
+      "real",
+      (MatrixXld(*)(const Eigen::MatrixBase<MatrixXcld> &))&real<MatrixXcld>);
+
+  bp::def("imag",
+          (MatrixXf(*)(const Eigen::MatrixBase<MatrixXcf> &))&imag<MatrixXcf>);
+  bp::def("imag",
+          (MatrixXd(*)(const Eigen::MatrixBase<MatrixXcd> &))&imag<MatrixXcd>);
+  bp::def(
+      "imag",
+      (MatrixXld(*)(const Eigen::MatrixBase<MatrixXcld> &))&imag<MatrixXcld>);
 }

unittest/deprecation_policy.cpp

+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/deprecation-policy.hpp"
+
+#include <iostream>
+
+namespace bp = boost::python;
+using eigenpy::DeprecationType;
+
+void some_deprecated_function() {
+  std::cout << "Calling this should produce a warning" << std::endl;
+}
+
+void some_future_deprecated_function() {
+  std::cout
+      << "Calling this should produce a warning about a future deprecation"
+      << std::endl;
+}
+
+class X {
+ public:
+  void deprecated_member_function() {}
+};
+
+BOOST_PYTHON_MODULE(deprecation_policy) {
+  bp::def("some_deprecated_function", some_deprecated_function,
+          eigenpy::deprecated_function<DeprecationType::DEPRECATION>());
+  bp::def("some_future_deprecated_function", some_future_deprecated_function,
+          eigenpy::deprecated_function<DeprecationType::FUTURE>());
+
+  bp::class_<X>("X", bp::init<>(bp::args("self")))
+      .def("deprecated_member_function", &X::deprecated_member_function,
+           eigenpy::deprecated_member<>());
+}

unittest/eigen_ref.cpp

 /*
  * Copyright 2014-2019, CNRS
- * Copyright 2018-2020, INRIA
+ * Copyright 2018-2022, INRIA
  */
 
-#include "eigenpy/eigenpy.hpp"
 #include <iostream>
 
+#include "eigenpy/eigenpy.hpp"
+#include "eigenpy/eigen-from-python.hpp"
+
 using namespace Eigen;
 using namespace eigenpy;
 
-template<typename MatType>
-void printMatrix(const Eigen::Ref<const MatType> mat)
-{
-  if(MatType::IsVectorAtCompileTime)
-    std::cout << "isVector" << std::endl;
-  std::cout << "size: cols " << mat.cols() << " rows " << mat.rows() << std::endl;
+template <typename MatType>
+void printMatrix(const Eigen::Ref<const MatType> mat) {
+  if (MatType::IsVectorAtCompileTime) std::cout << "isVector" << std::endl;
+  std::cout << "input size: cols " << mat.cols() << " rows " << mat.rows()
+            << std::endl;
   std::cout << mat << std::endl;
 }
 
-template<typename VecType>
-void printVector(const Eigen::Ref<const VecType> & vec)
-{
+template <typename VecType>
+void printVector(const Eigen::Ref<const VecType>& vec) {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VecType);
   printMatrix(vec);
 }
 
-template<typename MatType>
-void setOnes(Eigen::Ref<MatType> mat)
-{
+template <typename MatType>
+void setOnes(Eigen::Ref<MatType> mat) {
   mat.setOnes();
 }
 
-template<typename MatType>
-void fill(Eigen::Ref<MatType> mat, const typename MatType::Scalar & value)
-{
+template <typename VecType>
+VecType copyVectorFromConstRef(const Eigen::Ref<const VecType> vec) {
+  std::cout << "copyVectorFromConstRef::vec: " << vec.transpose() << std::endl;
+  return VecType(vec);
+}
+
+template <typename MatType>
+Eigen::Ref<MatType> getBlock(Eigen::Ref<MatType> mat, Eigen::DenseIndex i,
+                             Eigen::DenseIndex j, Eigen::DenseIndex n,
+                             Eigen::DenseIndex m) {
+  return mat.block(i, j, n, m);
+}
+
+template <typename MatType>
+Eigen::Ref<MatType> editBlock(Eigen::Ref<MatType> mat, Eigen::DenseIndex i,
+                              Eigen::DenseIndex j, Eigen::DenseIndex n,
+                              Eigen::DenseIndex m) {
+  typename Eigen::Ref<MatType>::BlockXpr B = mat.block(i, j, n, m);
+  int k = 0;
+  for (int i = 0; i < B.rows(); ++i) {
+    for (int j = 0; j < B.cols(); ++j) {
+      B(i, j) = k++;
+    }
+  }
+  std::cout << "B:\n" << B << std::endl;
+  return mat;
+}
+
+template <typename MatType>
+void fill(Eigen::Ref<MatType> mat, const typename MatType::Scalar& value) {
   mat.fill(value);
 }
 
-template<typename MatType>
-Eigen::Ref<MatType> asRef(const int rows, const int cols)
-{
-  static MatType mat(rows,cols);
-  std::cout << "mat:\n" << mat << std::endl;
+/// Get ref to a static matrix of size ( @p rows, @p cols )
+template <typename MatType>
+Eigen::Ref<MatType> getRefToStatic(const int rows, const int cols) {
+  static MatType mat(rows, cols);
+  std::cout << "create ref to matrix of size (" << rows << "," << cols << ")\n";
   return mat;
 }
 
-BOOST_PYTHON_MODULE(eigen_ref)
-{
+template <typename MatType>
+Eigen::Ref<MatType> asRef(Eigen::Ref<MatType> mat) {
+  std::cout << "create Ref to input mutable Ref\n";
+  return Eigen::Ref<MatType>(mat);
+}
+
+template <typename MatType>
+const Eigen::Ref<const MatType> asConstRef(Eigen::Ref<const MatType> mat) {
+  return Eigen::Ref<const MatType>(mat);
+}
+
+struct modify_block {
+  MatrixXd J;
+  modify_block() : J(10, 10) { J.setZero(); }
+  void modify(int n, int m) { call(J.topLeftCorner(n, m)); }
+  virtual void call(Eigen::Ref<MatrixXd> mat) = 0;
+};
+
+struct modify_block_wrap : modify_block, bp::wrapper<modify_block> {
+  modify_block_wrap() : modify_block() {}
+  void call(Eigen::Ref<MatrixXd> mat) { this->get_override("call")(mat); }
+};
+
+struct has_ref_member {
+  MatrixXd J;
+  Eigen::Ref<MatrixXd> Jref;
+  has_ref_member() : J(4, 4), Jref(J.topRightCorner(3, 3)) { J.setZero(); }
+};
+
+BOOST_PYTHON_MODULE(eigen_ref) {
   namespace bp = boost::python;
   eigenpy::enableEigenPy();
-  
+
   bp::def("printMatrix", printMatrix<Vector3d>);
   bp::def("printMatrix", printMatrix<VectorXd>);
   bp::def("printMatrix", printMatrix<MatrixXd>);
-  
+
   bp::def("printVector", printVector<VectorXd>);
   bp::def("printRowVector", printVector<RowVectorXd>);
 
   bp::def("setOnes", setOnes<Vector3d>);
   bp::def("setOnes", setOnes<VectorXd>);
   bp::def("setOnes", setOnes<MatrixXd>);
-  
+
   bp::def("fillVec3", fill<Vector3d>);
   bp::def("fillVec", fill<VectorXd>);
   bp::def("fill", fill<MatrixXd>);
-  
+
+  bp::def("getRefToStatic", getRefToStatic<MatrixXd>);
   bp::def("asRef", asRef<MatrixXd>);
+  bp::def("asConstRef", asConstRef<MatrixXd>);
+
+  bp::def("getBlock", &getBlock<MatrixXd>);
+  bp::def("editBlock", &editBlock<MatrixXd>);
+
+  bp::def("copyVectorFromConstRef", &copyVectorFromConstRef<VectorXd>);
+  bp::def("copyRowVectorFromConstRef", &copyVectorFromConstRef<RowVectorXd>);
+
+  bp::class_<modify_block_wrap, boost::noncopyable>("modify_block",
+                                                    bp::init<>())
+      .def_readonly("J", &modify_block::J)
+      .def("modify", &modify_block::modify)
+      .def("call", bp::pure_virtual(&modify_block_wrap::call));
+
+  bp::class_<has_ref_member, boost::noncopyable>("has_ref_member", bp::init<>())
+      .def_readonly("J", &has_ref_member::J)
+      .add_property(
+          "Jref",
+          bp::make_getter(&has_ref_member::Jref,
+                          bp::return_value_policy<bp::return_by_value>()));
+  // can't return Eigen::Ref by reference but by value
+  // (def_readonly creates a by-reference getter)
 }

unittest/geometry.cpp

 /*
  * Copyright 2014-2019, CNRS
- * Copyright 2018-2019, INRIA
+ * Copyright 2018-2023, INRIA
  */
 
 #include "eigenpy/eigenpy.hpp"
 #include "eigenpy/geometry.hpp"
-#include <Eigen/Geometry>
-#include <iostream>
- 
+
 namespace bp = boost::python;
 
-Eigen::AngleAxisd testOutAngleAxis()
-{
-  return Eigen::AngleAxisd(.1,Eigen::Vector3d::UnitZ());
+Eigen::AngleAxisd testOutAngleAxis() {
+  return Eigen::AngleAxisd(.1, Eigen::Vector3d::UnitZ());
 }
 
-double testInAngleAxis(Eigen::AngleAxisd aa)
-{
-  return aa.angle();
-}
+double testInAngleAxis(Eigen::AngleAxisd aa) { return aa.angle(); }
 
-Eigen::Quaterniond testOutQuaternion()
-{
-  Eigen::Quaterniond res(1,2,3,4);
+Eigen::Quaterniond testOutQuaternion() {
+  Eigen::Quaterniond res(1, 2, 3, 4);
   return res;
 }
-double testInQuaternion( Eigen::Quaterniond q )
-{
-  return q.norm(); 
-}
+double testInQuaternion(Eigen::Quaterniond q) { return q.norm(); }
 
-BOOST_PYTHON_MODULE(geometry)
-{
+BOOST_PYTHON_MODULE(geometry) {
   eigenpy::enableEigenPy();
 
   eigenpy::exposeAngleAxis();
   eigenpy::exposeQuaternion();
 
-  bp::def("testOutAngleAxis",&testOutAngleAxis);
-  bp::def("testInAngleAxis",&testInAngleAxis);
-
-  bp::def("testOutQuaternion",&testOutQuaternion);
-  bp::def("testInQuaternion",&testInQuaternion);
+  bp::def("testOutAngleAxis", &testOutAngleAxis);
+  bp::def("testInAngleAxis", &testInAngleAxis);
 
+  bp::def("testOutQuaternion", &testOutQuaternion);
+  bp::def("testInQuaternion", &testInQuaternion);
 }
-

unittest/matrix.cpp

 /*
- * Copyright 2014-2019, CNRS
- * Copyright 2018-2020, INRIA
+ * Copyright 2014-2022 CNRS INRIA
  */
 
-#include "eigenpy/eigenpy.hpp"
 #include <iostream>
 
-template<typename Scalar>
-Eigen::Matrix<Scalar,Eigen::Dynamic,1> vector1x1(const Scalar & value)
-{
-  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> ReturnType;
-  return ReturnType::Constant(1,value);
-} 
+#include "eigenpy/eigenpy.hpp"
+
+template <typename Scalar>
+Eigen::Matrix<Scalar, Eigen::Dynamic, 1> vector1x1(const Scalar& value) {
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> ReturnType;
+  return ReturnType::Constant(1, value);
+}
 
-template<typename Scalar>
-Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> matrix1x1(const Scalar & value)
-{
-  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> ReturnType; 
-  return ReturnType::Constant(1,1,value);
+template <typename Scalar>
+Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> matrix1x1(
+    const Scalar& value) {
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> ReturnType;
+  return ReturnType::Constant(1, 1, value);
 }
 
-Eigen::VectorXd emptyVector()
-{
+template <typename Scalar>
+void matrix1x1_input(const Eigen::Matrix<Scalar, 1, 1>& mat) {
+  std::cout << mat << std::endl;
+}
+
+Eigen::VectorXd emptyVector() {
   Eigen::VectorXd vec;
   vec.resize(0);
   return vec;
 }
 
-Eigen::MatrixXd emptyMatrix()
-{
-  return Eigen::MatrixXd(0,0);
-}
+Eigen::MatrixXd emptyMatrix() { return Eigen::MatrixXd(0, 0); }
 
-Eigen::MatrixXd naturals(int R,int C,bool verbose)
-{
-  Eigen::MatrixXd mat(R,C);
-  for(int r=0;r<R;++r)
-    for(int c=0;c<C;++c)
-      mat(r,c) = r*C+c;
+Eigen::MatrixXd naturals(int R, int C, bool verbose) {
+  Eigen::MatrixXd mat(R, C);
+  for (int r = 0; r < R; ++r)
+    for (int c = 0; c < C; ++c) mat(r, c) = r * C + c;
 
-  if(verbose)
-    std::cout << "EigenMat = " << mat << std::endl;
+  if (verbose) std::cout << "EigenMat = " << mat << std::endl;
   return mat;
 }
 
-Eigen::VectorXd naturals(int R,bool verbose)
-{
+Eigen::VectorXd naturals(int R, bool verbose) {
   Eigen::VectorXd mat(R);
-  for(int r=0;r<R;++r)
-    mat[r] = r;
+  for (int r = 0; r < R; ++r) mat[r] = r;
 
-  if(verbose)
-    std::cout << "EigenMat = " << mat << std::endl;
+  if (verbose) std::cout << "EigenMat = " << mat << std::endl;
   return mat;
 }
 
-Eigen::Matrix3d naturals(bool verbose)
-{
+Eigen::Matrix3d naturals(bool verbose) {
   Eigen::Matrix3d mat;
-  for(int r=0;r<3;++r)
-    for(int c=0;c<3;++c)
-      mat(r,c) = r*3+c;
+  for (int r = 0; r < 3; ++r)
+    for (int c = 0; c < 3; ++c) mat(r, c) = r * 3 + c;
 
-  if(verbose)
-    std::cout << "EigenMat = " << mat << std::endl;
+  if (verbose) std::cout << "EigenMat = " << mat << std::endl;
   return mat;
 }
 
-template<typename MatType>
-Eigen::MatrixXd reflex(const MatType & M, bool verbose)
-{
-  if(verbose)
-    std::cout << "EigenMat = " << M << std::endl;
+template <typename MatType>
+Eigen::MatrixXd reflex(const MatType& M, bool verbose) {
+  if (verbose) std::cout << "EigenMat = " << M << std::endl;
   return Eigen::MatrixXd(M);
 }
 
-template<typename MatrixDerived>
-MatrixDerived base(const Eigen::MatrixBase<MatrixDerived> & m)
-{
+template <typename MatrixDerived>
+MatrixDerived base(const Eigen::MatrixBase<MatrixDerived>& m) {
   return m.derived();
 }
 
-template<typename MatrixDerived>
-MatrixDerived plain(const Eigen::PlainObjectBase<MatrixDerived> & m)
-{
+template <typename MatrixDerived>
+MatrixDerived plain(const Eigen::PlainObjectBase<MatrixDerived>& m) {
   return m.derived();
 }
 
-template<typename Scalar>
-Eigen::Matrix<Scalar,6,6> matrix6(const Scalar & value)
-{
-  typedef Eigen::Matrix<Scalar,6,6> ReturnType;
+template <typename Scalar>
+Eigen::Matrix<Scalar, 6, 6> matrix6(const Scalar& value) {
+  typedef Eigen::Matrix<Scalar, 6, 6> ReturnType;
   return ReturnType::Constant(value);
 }
 
-BOOST_PYTHON_MODULE(matrix)
-{
+template <typename Scalar>
+Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+generateRowMajorMatrix(const Eigen::DenseIndex rows,
+                       const Eigen::DenseIndex cols) {
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+      RowMajorMatrix;
+  RowMajorMatrix A(rows, cols);
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
+  Eigen::Map<Vector>(A.data(), A.size()) =
+      Vector::LinSpaced(A.size(), 1, static_cast<Scalar>(A.size()));
+  std::cout << "Matrix values:\n" << A << std::endl;
+  return A;
+}
+
+template <typename Scalar>
+Eigen::Matrix<Scalar, 1, Eigen::Dynamic, Eigen::RowMajor>
+generateRowMajorVector(const Eigen::DenseIndex size) {
+  typedef Eigen::Matrix<Scalar, 1, Eigen::Dynamic, Eigen::RowMajor>
+      RowMajorVector;
+  RowMajorVector A(size);
+  A.setLinSpaced(size, 1, static_cast<Scalar>(size));
+  std::cout << "Vector values: " << A.transpose() << std::endl;
+  return A;
+}
+
+template <typename Scalar>
+Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> generateColMajorMatrix(
+    const Eigen::DenseIndex rows, const Eigen::DenseIndex cols) {
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> ColMajorMatrix;
+  ColMajorMatrix A(rows, cols);
+  typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
+  Eigen::Map<Vector>(A.data(), A.size()) =
+      Vector::LinSpaced(A.size(), 1, static_cast<Scalar>(A.size()));
+  std::cout << "Matrix values:\n" << A << std::endl;
+  return A;
+}
+
+template <typename Scalar>
+Eigen::Matrix<Scalar, 1, Eigen::Dynamic> generateColMajorVector(
+    const Eigen::DenseIndex size) {
+  typedef Eigen::Matrix<Scalar, 1, Eigen::Dynamic> ColMajorVector;
+  ColMajorVector A(size);
+  A.setLinSpaced(size, 1, static_cast<Scalar>(size));
+  std::cout << "Vector values: " << A.transpose() << std::endl;
+  return A;
+}
+
+template <typename Matrix, typename ReturnMatrix>
+ReturnMatrix copy(const Eigen::MatrixBase<Matrix>& mat) {
+  return mat;
+}
+
+template <typename Matrix>
+Matrix copy_same(const Eigen::MatrixBase<Matrix>& mat) {
+  return mat;
+}
+
+BOOST_PYTHON_MODULE(matrix) {
   using namespace Eigen;
   namespace bp = boost::python;
   eigenpy::enableEigenPy();
 
   // Square matrix
-  typedef Eigen::Matrix<double,6,6> Matrix6;
+  typedef Eigen::Matrix<double, 6, 6> Matrix6;
   eigenpy::enableEigenPySpecific<Matrix6>();
-  
+
   // Non-square matrix
-  typedef Eigen::Matrix<double,4,6> Matrix46;
+  typedef Eigen::Matrix<double, 4, 6> Matrix46;
   eigenpy::enableEigenPySpecific<Matrix46>();
 
-  Eigen::MatrixXd (*naturalsXX)(int,int,bool) = naturals;
-  Eigen::VectorXd (*naturalsX)(int,bool) = naturals;
+  Eigen::MatrixXd (*naturalsXX)(int, int, bool) = naturals;
+  Eigen::VectorXd (*naturalsX)(int, bool) = naturals;
   Eigen::Matrix3d (*naturals33)(bool) = naturals;
 
   bp::def("vector1x1", vector1x1<double>);
   bp::def("matrix1x1", matrix1x1<double>);
+  bp::def("matrix1x1", matrix1x1_input<double>);
+  bp::def("matrix1x1_int", matrix1x1_input<int>);
 
   bp::def("naturals", naturalsXX);
   bp::def("naturalsX", naturalsX);
@@ -126,12 +170,64 @@ BOOST_PYTHON_MODULE(matrix)
 
   bp::def("emptyVector", emptyVector);
   bp::def("emptyMatrix", emptyMatrix);
-  
+
   bp::def("base", base<VectorXd>);
   bp::def("base", base<MatrixXd>);
-  
+
   bp::def("plain", plain<VectorXd>);
   bp::def("plain", plain<MatrixXd>);
 
   bp::def("matrix6", matrix6<double>);
+
+  bp::def("generateRowMajorMatrix", generateRowMajorMatrix<double>);
+  bp::def("generateRowMajorVector", generateRowMajorVector<double>);
+
+  bp::def("generateColMajorMatrix", generateColMajorMatrix<double>);
+  bp::def("generateColMajorVector", generateColMajorVector<double>);
+
+  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
+      RowMajorMatrixXd;
+  bp::def("asRowMajorFromColMajorMatrix",
+          copy<Eigen::MatrixXd, RowMajorMatrixXd>);
+  bp::def("asRowMajorFromColMajorVector",
+          copy<Eigen::VectorXd, Eigen::RowVectorXd>);
+  bp::def("asRowMajorFromRowMajorMatrix",
+          copy<RowMajorMatrixXd, RowMajorMatrixXd>);
+  bp::def("asRowMajorFromRowMajorVector",
+          copy<Eigen::RowVectorXd, Eigen::RowVectorXd>);
+
+  bp::def("copyBoolToBool", copy_same<Eigen::Matrix<bool, -1, -1> >);
+
+  bp::def("copyInt8ToInt8", copy_same<Eigen::Matrix<int8_t, -1, -1> >);
+  bp::def("copyCharToChar", copy_same<Eigen::Matrix<char, -1, -1> >);
+  bp::def("copyUCharToUChar", copy_same<Eigen::Matrix<unsigned char, -1, -1> >);
+
+  bp::def("copyInt16ToInt16", copy_same<Eigen::Matrix<int16_t, -1, -1> >);
+  bp::def("copyUInt16ToUInt16", copy_same<Eigen::Matrix<uint16_t, -1, -1> >);
+
+  bp::def("copyInt32ToInt32", copy_same<Eigen::Matrix<int32_t, -1, -1> >);
+  bp::def("copyUInt32ToUInt32", copy_same<Eigen::Matrix<uint32_t, -1, -1> >);
+
+  bp::def("copyInt64ToInt64", copy_same<Eigen::Matrix<int64_t, -1, -1> >);
+  bp::def("copyUInt64ToUInt64", copy_same<Eigen::Matrix<uint64_t, -1, -1> >);
+
+  bp::def("copyLongToLong", copy_same<Eigen::Matrix<long, -1, -1> >);
+  bp::def("copyULongToULong", copy_same<Eigen::Matrix<unsigned long, -1, -1> >);
+
+  bp::def("copyLongLongToLongLong",
+          copy_same<Eigen::Matrix<long long, -1, -1> >);
+  bp::def("copyULongLongToULongLong",
+          copy_same<Eigen::Matrix<unsigned long long, -1, -1> >);
+
+  bp::def("copyFloatToFloat", copy_same<Eigen::Matrix<float, -1, -1> >);
+  bp::def("copyDoubleToDouble", copy_same<Eigen::Matrix<double, -1, -1> >);
+  bp::def("copyLongDoubleToLongDouble",
+          copy_same<Eigen::Matrix<long double, -1, -1> >);
+
+  bp::def("copyCFloatToCFloat",
+          copy_same<Eigen::Matrix<std::complex<float>, -1, -1> >);
+  bp::def("copyCDoubleToCDouble",
+          copy_same<Eigen::Matrix<std::complex<double>, -1, -1> >);
+  bp::def("copyCLongDoubleToCLongDouble",
+          copy_same<Eigen::Matrix<std::complex<long double>, -1, -1> >);
 }

unittest/multiple_registration.cpp

+#include "eigenpy/registration.hpp"
+#include <cstdio>
+
+namespace bp = boost::python;
+
+class X {
+ public:
+  X() {}
+  void operator()() { printf("DOOT\n"); }
+};
+
+class X_wrapper : public X, bp::wrapper<X> {
+ public:
+  static void expose() {
+    if (!eigenpy::register_symbolic_link_to_registered_type<X>()) {
+      bp::class_<X>("X", bp::init<>()).def("__call__", &X::operator());
+    }
+  }
+};
+
+BOOST_PYTHON_MODULE(multiple_registration) {
+  X_wrapper::expose();
+  X_wrapper::expose();
+  X_wrapper::expose();
+  X_wrapper::expose();
+  X_wrapper::expose();
+  X_wrapper::expose();
+}

unittest/packaging/cmake-2x/extra_lib.cpp

 
 #include <eigenpy/version.hpp>
 
-int main(int /*argc*/, char ** /*argv*/)
-{
-  eigenpy::checkVersionAtLeast(0,0,0);
+int main(int /*argc*/, char** /*argv*/) {
+  eigenpy::checkVersionAtLeast(0, 0, 0);
   return 0;
 }

unittest/packaging/cmake/CMakeLists.txt

@@ -2,8 +2,6 @@ cmake_minimum_required(VERSION 2.6)
 
 project(ExtraLib CXX)
 find_package(eigenpy REQUIRED)
-find_package(PythonInterp REQUIRED)
-find_package(PythonLibs REQUIRED)
 
 add_executable(extra_lib extra_lib.cpp)
-target_link_libraries(extra_lib PUBLIC eigenpy::eigenpy ${PYTHON_LIBRARIES})
+target_link_libraries(extra_lib PUBLIC eigenpy::eigenpy)

unittest/packaging/cmake/extra_lib.cpp

 
 #include <eigenpy/version.hpp>
 
-int main(int /*argc*/, char ** /*argv*/)
-{
-  eigenpy::checkVersionAtLeast(0,0,0);
+int main(int /*argc*/, char** /*argv*/) {
+  eigenpy::checkVersionAtLeast(0, 0, 0);
   return 0;
 }

unittest/packaging/pkgconfig/extra_lib.cpp

 
 #include <eigenpy/version.hpp>
 
-int main(int /*argc*/, char ** /*argv*/)
-{
-  eigenpy::checkVersionAtLeast(0,0,0);
+int main(int /*argc*/, char** /*argv*/) {
+  eigenpy::checkVersionAtLeast(0, 0, 0);
   return 0;
 }

unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLDLT.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSimplicialLDLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)

unittest/python/decompositions/sparse/cholmod/test_CholmodSimplicialLLT.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSimplicialLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)

unittest/python/decompositions/sparse/cholmod/test_CholmodSupernodalLLT.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.CholmodSupernodalLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)

unittest/python/decompositions/sparse/test_Accelerate.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+rng = np.random.default_rng()
+
+
+def test(SolverType: type):
+    dim = 100
+    A = rng.random((dim, dim))
+    A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+    A = csc_matrix(A)
+
+    llt = SolverType(A)
+
+    assert llt.info() == eigenpy.ComputationInfo.Success
+
+    X = rng.random((dim, 20))
+    B = A.dot(X)
+    X_est = llt.solve(B)
+    #    import pdb; pdb.set_trace()
+    assert eigenpy.is_approx(X, X_est)
+    assert eigenpy.is_approx(A.dot(X_est), B)
+
+    llt.analyzePattern(A)
+    llt.factorize(A)
+
+
+test(eigenpy.AccelerateLLT)
+test(eigenpy.AccelerateLDLT)
+test(eigenpy.AccelerateLDLTUnpivoted)
+test(eigenpy.AccelerateLDLTSBK)
+test(eigenpy.AccelerateLDLTTPP)
+test(eigenpy.AccelerateQR)
+# test(eigenpy.AccelerateCholeskyAtA) # This test is not passing. Seems there is a bug in Eigen with the support of Accelerate.

unittest/python/decompositions/sparse/test_SimplicialLDLT.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+A = csc_matrix(A)
+
+ldlt = eigenpy.SimplicialLDLT(A)
+
+assert ldlt.info() == eigenpy.ComputationInfo.Success
+
+L = ldlt.matrixL()
+U = ldlt.matrixU()
+D = csc_matrix(np.diag(ldlt.vectorD()))
+
+LDU = L @ D @ U
+assert eigenpy.is_approx(LDU.toarray(), A.toarray())
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = ldlt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+ldlt.analyzePattern(A)
+ldlt.factorize(A)
+permutation = ldlt.permutationP()

unittest/python/decompositions/sparse/test_SimplicialLLT.py

+import numpy as np
+from scipy.sparse import csc_matrix
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
+
+A = csc_matrix(A)
+
+llt = eigenpy.SimplicialLLT(A)
+
+assert llt.info() == eigenpy.ComputationInfo.Success
+
+L = llt.matrixL()
+U = llt.matrixU()
+
+LU = L @ U
+assert eigenpy.is_approx(LU.toarray(), A.toarray())
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)
+
+llt.analyzePattern(A)
+llt.factorize(A)
+permutation = llt.permutationP()

unittest/python/test_LDLT.py

-import eigenpy
-eigenpy.switchToNumpyArray()
-
 import numpy as np
-import numpy.linalg as la
+
+import eigenpy
 
 dim = 100
-A = np.random.rand(dim,dim)
+rng = np.random.default_rng()
 
-A = (A + A.T)*0.5 + np.diag(10. + np.random.rand(dim))
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
 
 ldlt = eigenpy.LDLT(A)
 
-L = ldlt.matrixL() 
-D = ldlt.vectorD() 
-P = ldlt.transpositionsP() 
+L = ldlt.matrixL()
+D = ldlt.vectorD()
+P = ldlt.transpositionsP()
+
+assert eigenpy.is_approx(
+    np.transpose(P).dot(L.dot(np.diag(D).dot(np.transpose(L).dot(P)))), A
+)
 
-assert eigenpy.is_approx(np.transpose(P).dot(L.dot(np.diag(D).dot(np.transpose(L).dot(P)))),A)
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = ldlt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)

unittest/python/test_LLT.py

-import eigenpy
-eigenpy.switchToNumpyArray()
-
 import numpy as np
-import numpy.linalg as la
+
+import eigenpy
 
 dim = 100
-A = np.random.rand(dim,dim)
+rng = np.random.default_rng()
 
-A = (A + A.T)*0.5 + np.diag(10. + np.random.rand(dim))
+A = rng.random((dim, dim))
+A = (A + A.T) * 0.5 + np.diag(10.0 + rng.random(dim))
 
 llt = eigenpy.LLT(A)
 
-L = llt.matrixL() 
+L = llt.matrixL()
+assert eigenpy.is_approx(L.dot(np.transpose(L)), A)
 
-assert eigenpy.is_approx(L.dot(np.transpose(L)),A)
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = llt.solve(B)
+assert eigenpy.is_approx(X, X_est)
+assert eigenpy.is_approx(A.dot(X_est), B)

unittest/python/test_MINRES.py

+import numpy as np
+
+import eigenpy
+
+dim = 100
+rng = np.random.default_rng()
+A = np.eye(dim)
+
+minres = eigenpy.MINRES(A)
+
+X = rng.random((dim, 20))
+B = A.dot(X)
+X_est = minres.solve(B)
+print("A.dot(X_est):", A.dot(X_est))
+print("B:", B)
+assert eigenpy.is_approx(A.dot(X_est), B, 1e-6)