diff options
Diffstat (limited to 'eigen/doc/examples')
19 files changed, 438 insertions, 150 deletions
diff --git a/eigen/doc/examples/CMakeLists.txt b/eigen/doc/examples/CMakeLists.txt index 08cf8ef..f7a1905 100644 --- a/eigen/doc/examples/CMakeLists.txt +++ b/eigen/doc/examples/CMakeLists.txt @@ -14,3 +14,8 @@ foreach(example_src ${examples_SRCS}) ) 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/eigen/doc/examples/CustomizingEigen_Inheritance.cpp b/eigen/doc/examples/CustomizingEigen_Inheritance.cpp new file mode 100644 index 0000000..48df64e --- /dev/null +++ b/eigen/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/eigen/doc/examples/Cwise_erf.cpp b/eigen/doc/examples/Cwise_erf.cpp new file mode 100644 index 0000000..e7cd2c1 --- /dev/null +++ b/eigen/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/eigen/doc/examples/Cwise_erfc.cpp b/eigen/doc/examples/Cwise_erfc.cpp new file mode 100644 index 0000000..d8bb04c --- /dev/null +++ b/eigen/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/eigen/doc/examples/Cwise_lgamma.cpp b/eigen/doc/examples/Cwise_lgamma.cpp new file mode 100644 index 0000000..f1c4f50 --- /dev/null +++ b/eigen/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; +}
\ No newline at end of file diff --git a/eigen/doc/examples/MatrixBase_cwise_const.cpp b/eigen/doc/examples/MatrixBase_cwise_const.cpp deleted file mode 100644 index 23700e0..0000000 --- a/eigen/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/eigen/doc/examples/TutorialInplaceLU.cpp b/eigen/doc/examples/TutorialInplaceLU.cpp new file mode 100644 index 0000000..cb9c59b --- /dev/null +++ b/eigen/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/eigen/doc/examples/TutorialLinAlgInverseDeterminant.cpp b/eigen/doc/examples/TutorialLinAlgInverseDeterminant.cpp index 43970ff..14dde5b 100644 --- a/eigen/doc/examples/TutorialLinAlgInverseDeterminant.cpp +++ b/eigen/doc/examples/TutorialLinAlgInverseDeterminant.cpp @@ -13,4 +13,4 @@ int main() 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/eigen/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_operatornorm.cpp b/eigen/doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_operatornorm.cpp new file mode 100644 index 0000000..62e28fc --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp b/eigen/doc/examples/make_circulant.cpp new file mode 100644 index 0000000..92e6aaa --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp.entry b/eigen/doc/examples/make_circulant.cpp.entry new file mode 100644 index 0000000..f9d2eb8 --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp.evaluator b/eigen/doc/examples/make_circulant.cpp.evaluator new file mode 100644 index 0000000..2ba79e7 --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp.expression b/eigen/doc/examples/make_circulant.cpp.expression new file mode 100644 index 0000000..380cd44 --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp.main b/eigen/doc/examples/make_circulant.cpp.main new file mode 100644 index 0000000..877f97f --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant.cpp.preamble b/eigen/doc/examples/make_circulant.cpp.preamble new file mode 100644 index 0000000..e575cce --- /dev/null +++ b/eigen/doc/examples/make_circulant.cpp.preamble @@ -0,0 +1,4 @@ +#include <Eigen/Core> +#include <iostream> + +template <class ArgType> class Circulant; diff --git a/eigen/doc/examples/make_circulant.cpp.traits b/eigen/doc/examples/make_circulant.cpp.traits new file mode 100644 index 0000000..4e04535 --- /dev/null +++ b/eigen/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/eigen/doc/examples/make_circulant2.cpp b/eigen/doc/examples/make_circulant2.cpp new file mode 100644 index 0000000..95d3dd3 --- /dev/null +++ b/eigen/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/eigen/doc/examples/matrixfree_cg.cpp b/eigen/doc/examples/matrixfree_cg.cpp index f0631c3..6a205ae 100644 --- a/eigen/doc/examples/matrixfree_cg.cpp +++ b/eigen/doc/examples/matrixfree_cg.cpp @@ -2,179 +2,127 @@ #include <Eigen/Core> #include <Eigen/Dense> #include <Eigen/IterativeLinearSolvers> +#include <unsupported/Eigen/IterativeSolvers> class MatrixReplacement; -template<typename Rhs> class MatrixReplacement_ProductReturnType; +using Eigen::SparseMatrix; namespace Eigen { namespace internal { + // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits: template<> - struct traits<MatrixReplacement> : Eigen::internal::traits<Eigen::SparseMatrix<double> > + struct traits<MatrixReplacement> : public Eigen::internal::traits<Eigen::SparseMatrix<double> > {}; - - template <typename Rhs> - struct traits<MatrixReplacement_ProductReturnType<Rhs> > { - // The equivalent plain objet type of the product. This type is used if the product needs to be evaluated into a temporary. - typedef Eigen::Matrix<typename Rhs::Scalar, Eigen::Dynamic, Rhs::ColsAtCompileTime> ReturnType; - }; } } -// Inheriting EigenBase should not be needed in the future. +// 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: - // Expose some compile-time information to Eigen: + // Required typedefs, constants, and method: typedef double Scalar; typedef double RealScalar; + typedef int StorageIndex; enum { ColsAtCompileTime = Eigen::Dynamic, - RowsAtCompileTime = Eigen::Dynamic, MaxColsAtCompileTime = Eigen::Dynamic, - MaxRowsAtCompileTime = Eigen::Dynamic + IsRowMajor = false }; - Index rows() const { return 4; } - Index cols() const { return 4; } + Index rows() const { return mp_mat->rows(); } + Index cols() const { return mp_mat->cols(); } - void resize(Index a_rows, Index a_cols) - { - // This method should not be needed in the future. - assert(a_rows==0 && a_cols==0 || a_rows==rows() && a_cols==cols()); - } - - // In the future, the return type should be Eigen::Product<MatrixReplacement,Rhs> template<typename Rhs> - MatrixReplacement_ProductReturnType<Rhs> operator*(const Eigen::MatrixBase<Rhs>& x) const { - return MatrixReplacement_ProductReturnType<Rhs>(*this, x.derived()); + 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) {} -// The proxy class representing the product of a MatrixReplacement with a MatrixBase<> -template<typename Rhs> -class MatrixReplacement_ProductReturnType : public Eigen::ReturnByValue<MatrixReplacement_ProductReturnType<Rhs> > { -public: - typedef MatrixReplacement::Index Index; - - // The ctor store references to the matrix and right-hand-side object (usually a vector). - MatrixReplacement_ProductReturnType(const MatrixReplacement& matrix, const Rhs& rhs) - : m_matrix(matrix), m_rhs(rhs) - {} - - Index rows() const { return m_matrix.rows(); } - Index cols() const { return m_rhs.cols(); } - - // This function is automatically called by Eigen. It must evaluate the product of matrix * rhs into y. - template<typename Dest> - void evalTo(Dest& y) const - { - y.setZero(4); - - y(0) += 2 * m_rhs(0); y(1) += 1 * m_rhs(0); - y(0) += 1 * m_rhs(1); y(1) += 2 * m_rhs(1); y(2) += 1 * m_rhs(1); - y(1) += 1 * m_rhs(2); y(2) += 2 * m_rhs(2); y(3) += 1 * m_rhs(2); - y(2) += 1 * m_rhs(3); y(3) += 2 * m_rhs(3); + void attachMyMatrix(const SparseMatrix<double> &mat) { + mp_mat = &mat; } + const SparseMatrix<double> my_matrix() const { return *mp_mat; } -protected: - const MatrixReplacement& m_matrix; - typename Rhs::Nested m_rhs; +private: + const SparseMatrix<double> *mp_mat; }; -/*****/ - -// This class simply warp a diagonal matrix as a Jacobi preconditioner. -// In the future such simple and generic wrapper should be shipped within Eigen itsel. -template <typename _Scalar> -class MyJacobiPreconditioner -{ - typedef _Scalar Scalar; - typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> Vector; - typedef typename Vector::Index Index; - - public: - // this typedef is only to export the scalar type and compile-time dimensions to solve_retval - typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixType; - - MyJacobiPreconditioner() : m_isInitialized(false) {} - - void setInvDiag(const Eigen::VectorXd &invdiag) { - m_invdiag=invdiag; - m_isInitialized=true; - } - - Index rows() const { return m_invdiag.size(); } - Index cols() const { return m_invdiag.size(); } - - template<typename MatType> - MyJacobiPreconditioner& analyzePattern(const MatType& ) { return *this; } - - template<typename MatType> - MyJacobiPreconditioner& factorize(const MatType& mat) { return *this; } - - template<typename MatType> - MyJacobiPreconditioner& compute(const MatType& mat) { return *this; } - - template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const - { - x = m_invdiag.array() * b.array() ; - } - - template<typename Rhs> inline const Eigen::internal::solve_retval<MyJacobiPreconditioner, Rhs> - solve(const Eigen::MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "MyJacobiPreconditioner is not initialized."); - eigen_assert(m_invdiag.size()==b.rows() - && "MyJacobiPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); - return Eigen::internal::solve_retval<MyJacobiPreconditioner, Rhs>(*this, b.derived()); - } - - protected: - Vector m_invdiag; - bool m_isInitialized; -}; - +// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl: namespace Eigen { namespace internal { -template<typename _MatrixType, typename Rhs> -struct solve_retval<MyJacobiPreconditioner<_MatrixType>, Rhs> - : solve_retval_base<MyJacobiPreconditioner<_MatrixType>, Rhs> -{ - typedef MyJacobiPreconditioner<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const + 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> > { - dec()._solve(rhs(),dst); - } -}; + 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"); + + // 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; - Eigen::VectorXd b(4), x; - b << 1, 1, 1, 1; + A.attachMyMatrix(S); - // solve Ax = b using CG with matrix-free version: - Eigen::ConjugateGradient < MatrixReplacement, Eigen::Lower|Eigen::Upper, MyJacobiPreconditioner<double> > cg; + Eigen::VectorXd b(n), x; + b.setRandom(); - Eigen::VectorXd invdiag(4); - invdiag << 1./3., 1./4., 1./4., 1./3.; + // 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; + } - cg.preconditioner().setInvDiag(invdiag); - cg.compute(A); - x = cg.solve(b); + { + 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; + } - std::cout << "#iterations: " << cg.iterations() << std::endl; - std::cout << "estimated error: " << cg.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/eigen/doc/examples/nullary_indexing.cpp b/eigen/doc/examples/nullary_indexing.cpp new file mode 100644 index 0000000..e27c358 --- /dev/null +++ b/eigen/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"; +} + |