summaryrefslogtreecommitdiffhomepage
path: root/eigen/doc/examples/matrixfree_cg.cpp
diff options
context:
space:
mode:
authorStanislaw Halik <sthalik@misaki.pl>2019-03-03 21:09:10 +0100
committerStanislaw Halik <sthalik@misaki.pl>2019-03-03 21:10:13 +0100
commitf0238cfb6997c4acfc2bd200de7295f3fa36968f (patch)
treeb215183760e4f615b9c1dabc1f116383b72a1b55 /eigen/doc/examples/matrixfree_cg.cpp
parent543edd372a5193d04b3de9f23c176ab439e51b31 (diff)
don't index Eigen
Diffstat (limited to 'eigen/doc/examples/matrixfree_cg.cpp')
-rw-r--r--eigen/doc/examples/matrixfree_cg.cpp129
1 files changed, 0 insertions, 129 deletions
diff --git a/eigen/doc/examples/matrixfree_cg.cpp b/eigen/doc/examples/matrixfree_cg.cpp
deleted file mode 100644
index 7469938..0000000
--- a/eigen/doc/examples/matrixfree_cg.cpp
+++ /dev/null
@@ -1,129 +0,0 @@
-#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;
- }
-}