diff options
| author | Stanislaw Halik <sthalik@misaki.pl> | 2019-03-03 21:09:10 +0100 |
|---|---|---|
| committer | Stanislaw Halik <sthalik@misaki.pl> | 2019-03-03 21:10:13 +0100 |
| commit | f0238cfb6997c4acfc2bd200de7295f3fa36968f (patch) | |
| tree | b215183760e4f615b9c1dabc1f116383b72a1b55 /eigen/doc/examples/matrixfree_cg.cpp | |
| parent | 543edd372a5193d04b3de9f23c176ab439e51b31 (diff) | |
don't index Eigen
Diffstat (limited to 'eigen/doc/examples/matrixfree_cg.cpp')
| -rw-r--r-- | eigen/doc/examples/matrixfree_cg.cpp | 129 |
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; - } -} |
