Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/doc/examples/CMakeLists.txt b/doc/examples/CMakeLists.txt
index 08cf8ef..f7a1905 100644
--- a/doc/examples/CMakeLists.txt
+++ b/doc/examples/CMakeLists.txt
@@ -14,3 +14,8 @@
)
add_dependencies(all_examples ${example})
endforeach(example_src)
+
+check_cxx_compiler_flag("-std=c++11" EIGEN_COMPILER_SUPPORT_CPP11)
+if(EIGEN_COMPILER_SUPPORT_CPP11)
+ei_add_target_property(nullary_indexing COMPILE_FLAGS "-std=c++11")
+endif()
\ No newline at end of file
diff --git a/doc/examples/CustomizingEigen_Inheritance.cpp b/doc/examples/CustomizingEigen_Inheritance.cpp
new file mode 100644
index 0000000..48df64e
--- /dev/null
+++ b/doc/examples/CustomizingEigen_Inheritance.cpp
@@ -0,0 +1,30 @@
+#include <Eigen/Core>
+#include <iostream>
+
+class MyVectorType : public Eigen::VectorXd
+{
+public:
+ MyVectorType(void):Eigen::VectorXd() {}
+
+ // This constructor allows you to construct MyVectorType from Eigen expressions
+ template<typename OtherDerived>
+ MyVectorType(const Eigen::MatrixBase<OtherDerived>& other)
+ : Eigen::VectorXd(other)
+ { }
+
+ // This method allows you to assign Eigen expressions to MyVectorType
+ template<typename OtherDerived>
+ MyVectorType& operator=(const Eigen::MatrixBase <OtherDerived>& other)
+ {
+ this->Eigen::VectorXd::operator=(other);
+ return *this;
+ }
+};
+
+int main()
+{
+ MyVectorType v = MyVectorType::Ones(4);
+ v(2) += 10;
+ v = 2 * v;
+ std::cout << v.transpose() << std::endl;
+}
diff --git a/doc/examples/Cwise_erf.cpp b/doc/examples/Cwise_erf.cpp
new file mode 100644
index 0000000..e7cd2c1
--- /dev/null
+++ b/doc/examples/Cwise_erf.cpp
@@ -0,0 +1,9 @@
+#include <Eigen/Core>
+#include <unsupported/Eigen/SpecialFunctions>
+#include <iostream>
+using namespace Eigen;
+int main()
+{
+ Array4d v(-0.5,2,0,-7);
+ std::cout << v.erf() << std::endl;
+}
diff --git a/doc/examples/Cwise_erfc.cpp b/doc/examples/Cwise_erfc.cpp
new file mode 100644
index 0000000..d8bb04c
--- /dev/null
+++ b/doc/examples/Cwise_erfc.cpp
@@ -0,0 +1,9 @@
+#include <Eigen/Core>
+#include <unsupported/Eigen/SpecialFunctions>
+#include <iostream>
+using namespace Eigen;
+int main()
+{
+ Array4d v(-0.5,2,0,-7);
+ std::cout << v.erfc() << std::endl;
+}
diff --git a/doc/examples/Cwise_lgamma.cpp b/doc/examples/Cwise_lgamma.cpp
new file mode 100644
index 0000000..6bfaccb
--- /dev/null
+++ b/doc/examples/Cwise_lgamma.cpp
@@ -0,0 +1,9 @@
+#include <Eigen/Core>
+#include <unsupported/Eigen/SpecialFunctions>
+#include <iostream>
+using namespace Eigen;
+int main()
+{
+ Array4d v(0.5,10,0,-1);
+ std::cout << v.lgamma() << std::endl;
+}
diff --git a/doc/examples/MatrixBase_cwise_const.cpp b/doc/examples/MatrixBase_cwise_const.cpp
deleted file mode 100644
index 23700e0..0000000
--- a/doc/examples/MatrixBase_cwise_const.cpp
+++ /dev/null
@@ -1,18 +0,0 @@
-#define EIGEN2_SUPPORT
-#include <Eigen/Core>
-#include <iostream>
-
-using namespace Eigen;
-using namespace std;
-
-int main()
-{
- Matrix3i m = Matrix3i::Random();
- cout << "Here is the matrix m:" << endl << m << endl;
- Matrix3i n = Matrix3i::Random();
- cout << "And here is the matrix n:" << endl << n << endl;
- cout << "The coefficient-wise product of m and n is:" << endl;
- cout << m.cwise() * n << endl;
- cout << "Taking the cube of the coefficients of m yields:" << endl;
- cout << m.cwise().pow(3) << endl;
-}
diff --git a/doc/examples/TutorialInplaceLU.cpp b/doc/examples/TutorialInplaceLU.cpp
new file mode 100644
index 0000000..cb9c59b
--- /dev/null
+++ b/doc/examples/TutorialInplaceLU.cpp
@@ -0,0 +1,61 @@
+#include <iostream>
+struct init {
+ init() { std::cout << "[" << "init" << "]" << std::endl; }
+};
+init init_obj;
+// [init]
+#include <iostream>
+#include <Eigen/Dense>
+
+using namespace std;
+using namespace Eigen;
+
+int main()
+{
+ MatrixXd A(2,2);
+ A << 2, -1, 1, 3;
+ cout << "Here is the input matrix A before decomposition:\n" << A << endl;
+cout << "[init]" << endl;
+
+cout << "[declaration]" << endl;
+ PartialPivLU<Ref<MatrixXd> > lu(A);
+ cout << "Here is the input matrix A after decomposition:\n" << A << endl;
+cout << "[declaration]" << endl;
+
+cout << "[matrixLU]" << endl;
+ cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl;
+cout << "[matrixLU]" << endl;
+
+cout << "[solve]" << endl;
+ MatrixXd A0(2,2); A0 << 2, -1, 1, 3;
+ VectorXd b(2); b << 1, 2;
+ VectorXd x = lu.solve(b);
+ cout << "Residual: " << (A0 * x - b).norm() << endl;
+cout << "[solve]" << endl;
+
+cout << "[modifyA]" << endl;
+ A << 3, 4, -2, 1;
+ x = lu.solve(b);
+ cout << "Residual: " << (A0 * x - b).norm() << endl;
+cout << "[modifyA]" << endl;
+
+cout << "[recompute]" << endl;
+ A0 = A; // save A
+ lu.compute(A);
+ x = lu.solve(b);
+ cout << "Residual: " << (A0 * x - b).norm() << endl;
+cout << "[recompute]" << endl;
+
+cout << "[recompute_bis0]" << endl;
+ MatrixXd A1(2,2);
+ A1 << 5,-2,3,4;
+ lu.compute(A1);
+ cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl;
+cout << "[recompute_bis0]" << endl;
+
+cout << "[recompute_bis1]" << endl;
+ x = lu.solve(b);
+ cout << "Residual: " << (A1 * x - b).norm() << endl;
+cout << "[recompute_bis1]" << endl;
+
+}
diff --git a/doc/examples/TutorialLinAlgInverseDeterminant.cpp b/doc/examples/TutorialLinAlgInverseDeterminant.cpp
index 43970ff..14dde5b 100644
--- a/doc/examples/TutorialLinAlgInverseDeterminant.cpp
+++ b/doc/examples/TutorialLinAlgInverseDeterminant.cpp
@@ -13,4 +13,4 @@
cout << "Here is the matrix A:\n" << A << endl;
cout << "The determinant of A is " << A.determinant() << endl;
cout << "The inverse of A is:\n" << A.inverse() << endl;
-}
\ No newline at end of file
+}
diff --git a/doc/examples/TutorialLinAlgSVDSolve.cpp b/doc/examples/TutorialLinAlgSVDSolve.cpp
index 9fbc031..f109f04 100644
--- a/doc/examples/TutorialLinAlgSVDSolve.cpp
+++ b/doc/examples/TutorialLinAlgSVDSolve.cpp
@@ -11,5 +11,5 @@
VectorXf b = VectorXf::Random(3);
cout << "Here is the right hand side b:\n" << b << endl;
cout << "The least-squares solution is:\n"
- << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
+ << A.bdcSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
}
diff --git a/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_operatornorm.cpp b/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_operatornorm.cpp
new file mode 100644
index 0000000..62e28fc
--- /dev/null
+++ b/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_operatornorm.cpp
@@ -0,0 +1,18 @@
+#include <Eigen/Dense>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+ MatrixXf m(2,2);
+ m << 1,-2,
+ -3,4;
+
+ cout << "1-norm(m) = " << m.cwiseAbs().colwise().sum().maxCoeff()
+ << " == " << m.colwise().lpNorm<1>().maxCoeff() << endl;
+
+ cout << "infty-norm(m) = " << m.cwiseAbs().rowwise().sum().maxCoeff()
+ << " == " << m.rowwise().lpNorm<1>().maxCoeff() << endl;
+}
diff --git a/doc/examples/Tutorial_simple_example_dynamic_size.cpp b/doc/examples/Tutorial_simple_example_dynamic_size.cpp
index 0f0280e..defcb1e 100644
--- a/doc/examples/Tutorial_simple_example_dynamic_size.cpp
+++ b/doc/examples/Tutorial_simple_example_dynamic_size.cpp
@@ -10,7 +10,7 @@
MatrixXi m(size,size+1); // a (size)x(size+1)-matrix of int's
for (int j=0; j<m.cols(); ++j) // loop over columns
for (int i=0; i<m.rows(); ++i) // loop over rows
- m(i,j) = i+j*m.rows(); // to access matrix coefficients,
+ m(i,j) = i+j*size; // to access matrix coefficients,
// use operator()(int,int)
std::cout << m << "\n\n";
}
diff --git a/doc/examples/make_circulant.cpp b/doc/examples/make_circulant.cpp
new file mode 100644
index 0000000..92e6aaa
--- /dev/null
+++ b/doc/examples/make_circulant.cpp
@@ -0,0 +1,11 @@
+/*
+This program is presented in several fragments in the doc page.
+Every fragment is in its own file; this file simply combines them.
+*/
+
+#include "make_circulant.cpp.preamble"
+#include "make_circulant.cpp.traits"
+#include "make_circulant.cpp.expression"
+#include "make_circulant.cpp.evaluator"
+#include "make_circulant.cpp.entry"
+#include "make_circulant.cpp.main"
diff --git a/doc/examples/make_circulant.cpp.entry b/doc/examples/make_circulant.cpp.entry
new file mode 100644
index 0000000..f9d2eb8
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.entry
@@ -0,0 +1,5 @@
+template <class ArgType>
+Circulant<ArgType> makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
+{
+ return Circulant<ArgType>(arg.derived());
+}
diff --git a/doc/examples/make_circulant.cpp.evaluator b/doc/examples/make_circulant.cpp.evaluator
new file mode 100644
index 0000000..2ba79e7
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.evaluator
@@ -0,0 +1,32 @@
+namespace Eigen {
+ namespace internal {
+ template<typename ArgType>
+ struct evaluator<Circulant<ArgType> >
+ : evaluator_base<Circulant<ArgType> >
+ {
+ typedef Circulant<ArgType> XprType;
+ typedef typename nested_eval<ArgType, XprType::ColsAtCompileTime>::type ArgTypeNested;
+ typedef typename remove_all<ArgTypeNested>::type ArgTypeNestedCleaned;
+ typedef typename XprType::CoeffReturnType CoeffReturnType;
+
+ enum {
+ CoeffReadCost = evaluator<ArgTypeNestedCleaned>::CoeffReadCost,
+ Flags = Eigen::ColMajor
+ };
+
+ evaluator(const XprType& xpr)
+ : m_argImpl(xpr.m_arg), m_rows(xpr.rows())
+ { }
+
+ CoeffReturnType coeff(Index row, Index col) const
+ {
+ Index index = row - col;
+ if (index < 0) index += m_rows;
+ return m_argImpl.coeff(index);
+ }
+
+ evaluator<ArgTypeNestedCleaned> m_argImpl;
+ const Index m_rows;
+ };
+ }
+}
diff --git a/doc/examples/make_circulant.cpp.expression b/doc/examples/make_circulant.cpp.expression
new file mode 100644
index 0000000..380cd44
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.expression
@@ -0,0 +1,20 @@
+template <class ArgType>
+class Circulant : public Eigen::MatrixBase<Circulant<ArgType> >
+{
+public:
+ Circulant(const ArgType& arg)
+ : m_arg(arg)
+ {
+ EIGEN_STATIC_ASSERT(ArgType::ColsAtCompileTime == 1,
+ YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
+ }
+
+ typedef typename Eigen::internal::ref_selector<Circulant>::type Nested;
+
+ typedef Eigen::Index Index;
+ Index rows() const { return m_arg.rows(); }
+ Index cols() const { return m_arg.rows(); }
+
+ typedef typename Eigen::internal::ref_selector<ArgType>::type ArgTypeNested;
+ ArgTypeNested m_arg;
+};
diff --git a/doc/examples/make_circulant.cpp.main b/doc/examples/make_circulant.cpp.main
new file mode 100644
index 0000000..877f97f
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.main
@@ -0,0 +1,8 @@
+int main()
+{
+ Eigen::VectorXd vec(4);
+ vec << 1, 2, 4, 8;
+ Eigen::MatrixXd mat;
+ mat = makeCirculant(vec);
+ std::cout << mat << std::endl;
+}
diff --git a/doc/examples/make_circulant.cpp.preamble b/doc/examples/make_circulant.cpp.preamble
new file mode 100644
index 0000000..e575cce
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.preamble
@@ -0,0 +1,4 @@
+#include <Eigen/Core>
+#include <iostream>
+
+template <class ArgType> class Circulant;
diff --git a/doc/examples/make_circulant.cpp.traits b/doc/examples/make_circulant.cpp.traits
new file mode 100644
index 0000000..4e04535
--- /dev/null
+++ b/doc/examples/make_circulant.cpp.traits
@@ -0,0 +1,19 @@
+namespace Eigen {
+ namespace internal {
+ template <class ArgType>
+ struct traits<Circulant<ArgType> >
+ {
+ typedef Eigen::Dense StorageKind;
+ typedef Eigen::MatrixXpr XprKind;
+ typedef typename ArgType::StorageIndex StorageIndex;
+ typedef typename ArgType::Scalar Scalar;
+ enum {
+ Flags = Eigen::ColMajor,
+ RowsAtCompileTime = ArgType::RowsAtCompileTime,
+ ColsAtCompileTime = ArgType::RowsAtCompileTime,
+ MaxRowsAtCompileTime = ArgType::MaxRowsAtCompileTime,
+ MaxColsAtCompileTime = ArgType::MaxRowsAtCompileTime
+ };
+ };
+ }
+}
diff --git a/doc/examples/make_circulant2.cpp b/doc/examples/make_circulant2.cpp
new file mode 100644
index 0000000..95d3dd3
--- /dev/null
+++ b/doc/examples/make_circulant2.cpp
@@ -0,0 +1,52 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+
+// [circulant_func]
+template<class ArgType>
+class circulant_functor {
+ const ArgType &m_vec;
+public:
+ circulant_functor(const ArgType& arg) : m_vec(arg) {}
+
+ const typename ArgType::Scalar& operator() (Index row, Index col) const {
+ Index index = row - col;
+ if (index < 0) index += m_vec.size();
+ return m_vec(index);
+ }
+};
+// [circulant_func]
+
+// [square]
+template<class ArgType>
+struct circulant_helper {
+ typedef Matrix<typename ArgType::Scalar,
+ ArgType::SizeAtCompileTime,
+ ArgType::SizeAtCompileTime,
+ ColMajor,
+ ArgType::MaxSizeAtCompileTime,
+ ArgType::MaxSizeAtCompileTime> MatrixType;
+};
+// [square]
+
+// [makeCirculant]
+template <class ArgType>
+CwiseNullaryOp<circulant_functor<ArgType>, typename circulant_helper<ArgType>::MatrixType>
+makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
+{
+ typedef typename circulant_helper<ArgType>::MatrixType MatrixType;
+ return MatrixType::NullaryExpr(arg.size(), arg.size(), circulant_functor<ArgType>(arg.derived()));
+}
+// [makeCirculant]
+
+// [main]
+int main()
+{
+ Eigen::VectorXd vec(4);
+ vec << 1, 2, 4, 8;
+ Eigen::MatrixXd mat;
+ mat = makeCirculant(vec);
+ std::cout << mat << std::endl;
+}
+// [main]
diff --git a/doc/examples/matrixfree_cg.cpp b/doc/examples/matrixfree_cg.cpp
new file mode 100644
index 0000000..7469938
--- /dev/null
+++ b/doc/examples/matrixfree_cg.cpp
@@ -0,0 +1,129 @@
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/IterativeLinearSolvers>
+#include <unsupported/Eigen/IterativeSolvers>
+
+class MatrixReplacement;
+using Eigen::SparseMatrix;
+
+namespace Eigen {
+namespace internal {
+ // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits:
+ template<>
+ struct traits<MatrixReplacement> : public Eigen::internal::traits<Eigen::SparseMatrix<double> >
+ {};
+}
+}
+
+// Example of a matrix-free wrapper from a user type to Eigen's compatible type
+// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix.
+class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
+public:
+ // Required typedefs, constants, and method:
+ typedef double Scalar;
+ typedef double RealScalar;
+ typedef int StorageIndex;
+ enum {
+ ColsAtCompileTime = Eigen::Dynamic,
+ MaxColsAtCompileTime = Eigen::Dynamic,
+ IsRowMajor = false
+ };
+
+ Index rows() const { return mp_mat->rows(); }
+ Index cols() const { return mp_mat->cols(); }
+
+ template<typename Rhs>
+ Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
+ return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
+ }
+
+ // Custom API:
+ MatrixReplacement() : mp_mat(0) {}
+
+ void attachMyMatrix(const SparseMatrix<double> &mat) {
+ mp_mat = &mat;
+ }
+ const SparseMatrix<double> my_matrix() const { return *mp_mat; }
+
+private:
+ const SparseMatrix<double> *mp_mat;
+};
+
+
+// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
+namespace Eigen {
+namespace internal {
+
+ template<typename Rhs>
+ struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
+ : generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
+ {
+ typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
+
+ template<typename Dest>
+ static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
+ {
+ // This method should implement "dst += alpha * lhs * rhs" inplace,
+ // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
+ assert(alpha==Scalar(1) && "scaling is not implemented");
+ EIGEN_ONLY_USED_FOR_DEBUG(alpha);
+
+ // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
+ // but let's do something fancier (and less efficient):
+ for(Index i=0; i<lhs.cols(); ++i)
+ dst += rhs(i) * lhs.my_matrix().col(i);
+ }
+ };
+
+}
+}
+
+int main()
+{
+ int n = 10;
+ Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
+ S = S.transpose()*S;
+
+ MatrixReplacement A;
+ A.attachMyMatrix(S);
+
+ Eigen::VectorXd b(n), x;
+ b.setRandom();
+
+ // Solve Ax = b using various iterative solver with matrix-free version:
+ {
+ Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
+ cg.compute(A);
+ x = cg.solve(b);
+ std::cout << "CG: #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
+ }
+
+ {
+ Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg;
+ bicg.compute(A);
+ x = bicg.solve(b);
+ std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
+ }
+
+ {
+ Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
+ gmres.compute(A);
+ x = gmres.solve(b);
+ std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
+ }
+
+ {
+ Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
+ gmres.compute(A);
+ x = gmres.solve(b);
+ std::cout << "DGMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
+ }
+
+ {
+ Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres;
+ minres.compute(A);
+ x = minres.solve(b);
+ std::cout << "MINRES: #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
+ }
+}
diff --git a/doc/examples/nullary_indexing.cpp b/doc/examples/nullary_indexing.cpp
new file mode 100644
index 0000000..e27c358
--- /dev/null
+++ b/doc/examples/nullary_indexing.cpp
@@ -0,0 +1,66 @@
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+
+// [functor]
+template<class ArgType, class RowIndexType, class ColIndexType>
+class indexing_functor {
+ const ArgType &m_arg;
+ const RowIndexType &m_rowIndices;
+ const ColIndexType &m_colIndices;
+public:
+ typedef Matrix<typename ArgType::Scalar,
+ RowIndexType::SizeAtCompileTime,
+ ColIndexType::SizeAtCompileTime,
+ ArgType::Flags&RowMajorBit?RowMajor:ColMajor,
+ RowIndexType::MaxSizeAtCompileTime,
+ ColIndexType::MaxSizeAtCompileTime> MatrixType;
+
+ indexing_functor(const ArgType& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
+ : m_arg(arg), m_rowIndices(row_indices), m_colIndices(col_indices)
+ {}
+
+ const typename ArgType::Scalar& operator() (Index row, Index col) const {
+ return m_arg(m_rowIndices[row], m_colIndices[col]);
+ }
+};
+// [functor]
+
+// [function]
+template <class ArgType, class RowIndexType, class ColIndexType>
+CwiseNullaryOp<indexing_functor<ArgType,RowIndexType,ColIndexType>, typename indexing_functor<ArgType,RowIndexType,ColIndexType>::MatrixType>
+indexing(const Eigen::MatrixBase<ArgType>& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
+{
+ typedef indexing_functor<ArgType,RowIndexType,ColIndexType> Func;
+ typedef typename Func::MatrixType MatrixType;
+ return MatrixType::NullaryExpr(row_indices.size(), col_indices.size(), Func(arg.derived(), row_indices, col_indices));
+}
+// [function]
+
+
+int main()
+{
+ std::cout << "[main1]\n";
+ Eigen::MatrixXi A = Eigen::MatrixXi::Random(4,4);
+ Array3i ri(1,2,1);
+ ArrayXi ci(6); ci << 3,2,1,0,0,2;
+ Eigen::MatrixXi B = indexing(A, ri, ci);
+ std::cout << "A =" << std::endl;
+ std::cout << A << std::endl << std::endl;
+ std::cout << "A([" << ri.transpose() << "], [" << ci.transpose() << "]) =" << std::endl;
+ std::cout << B << std::endl;
+ std::cout << "[main1]\n";
+
+ std::cout << "[main2]\n";
+ B = indexing(A, ri+1, ci);
+ std::cout << "A(ri+1,ci) =" << std::endl;
+ std::cout << B << std::endl << std::endl;
+#if __cplusplus >= 201103L
+ B = indexing(A, ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3));
+ std::cout << "A(ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3)) =" << std::endl;
+ std::cout << B << std::endl << std::endl;
+#endif
+ std::cout << "[main2]\n";
+}
+