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/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;
+ }
+}