src/IterativeLinearSolvers/BiCGSTAB.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_BICGSTAB_H
0012 #define EIGEN_BICGSTAB_H
0013
0014 namespace RivetEigen {
0015
0016 namespace internal {
0017
0018 /** \internal Low-level bi conjugate gradient stabilized algorithm
0019   * \param mat The matrix A
0020   * \param rhs The right hand side vector b
0021   * \param x On input and initial solution, on output the computed solution.
0022   * \param precond A preconditioner being able to efficiently solve for an
0023   *                approximation of Ax=b (regardless of b)
0024   * \param iters On input the max number of iteration, on output the number of performed iterations.
0025   * \param tol_error On input the tolerance error, on output an estimation of the relative error.
0026   * \return false in the case of numerical issue, for example a break down of BiCGSTAB.
0027   */
0028 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
0029 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
0030               const Preconditioner& precond, Index& iters,
0031               typename Dest::RealScalar& tol_error)
0032 {
0033   using std::sqrt;
0034   using std::abs;
0035   typedef typename Dest::RealScalar RealScalar;
0036   typedef typename Dest::Scalar Scalar;
0037   typedef Matrix<Scalar,Dynamic,1> VectorType;
0038   RealScalar tol = tol_error;
0039   Index maxIters = iters;
0040
0041   Index n = mat.cols();
0042   VectorType r  = rhs - mat * x;
0043   VectorType r0 = r;
0044
0045   RealScalar r0_sqnorm = r0.squaredNorm();
0046   RealScalar rhs_sqnorm = rhs.squaredNorm();
0047   if(rhs_sqnorm == 0)
0048   {
0049     x.setZero();
0050     return true;
0051   }
0052   Scalar rho    = 1;
0053   Scalar alpha  = 1;
0054   Scalar w      = 1;
0055
0056   VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
0057   VectorType y(n),  z(n);
0058   VectorType kt(n), ks(n);
0059
0060   VectorType s(n), t(n);
0061
0062   RealScalar tol2 = tol*tol*rhs_sqnorm;
0063   RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
0064   Index i = 0;
0065   Index restarts = 0;
0066
0067   while ( r.squaredNorm() > tol2 && i<maxIters )
0068   {
0069     Scalar rho_old = rho;
0070
0071     rho = r0.dot(r);
0072     if (abs(rho) < eps2*r0_sqnorm)
0073     {
0074       // The new residual vector became too orthogonal to the arbitrarily chosen direction r0
0075       // Let's restart with a new r0:
0076       r  = rhs - mat * x;
0077       r0 = r;
0078       rho = r0_sqnorm = r.squaredNorm();
0079       if(restarts++ == 0)
0080         i = 0;
0081     }
0082     Scalar beta = (rho/rho_old) * (alpha / w);
0083     p = r + beta * (p - w * v);
0084
0085     y = precond.solve(p);
0086
0087     v.noalias() = mat * y;
0088
0089     alpha = rho / r0.dot(v);
0090     s = r - alpha * v;
0091
0092     z = precond.solve(s);
0093     t.noalias() = mat * z;
0094
0095     RealScalar tmp = t.squaredNorm();
0096     if(tmp>RealScalar(0))
0097       w = t.dot(s) / tmp;
0098     else
0099       w = Scalar(0);
0100     x += alpha * y + w * z;
0101     r = s - w * t;
0102     ++i;
0103   }
0104   tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
0105   iters = i;
0106   return true;
0107 }
0108
0109 }
0110
0111 template< typename _MatrixType,
0112           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
0113 class BiCGSTAB;
0114
0115 namespace internal {
0116
0117 template< typename _MatrixType, typename _Preconditioner>
0118 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
0119 {
0120   typedef _MatrixType MatrixType;
0121   typedef _Preconditioner Preconditioner;
0122 };
0123
0124 }
0125
0126 /** \ingroup IterativeLinearSolvers_Module
0127   * \brief A bi conjugate gradient stabilized solver for sparse square problems
0128   *
0129   * This class allows to solve for A.x = b sparse linear problems using a bi conjugate gradient
0130   * stabilized algorithm. The vectors x and b can be either dense or sparse.
0131   *
0132   * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
0133   * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
0134   *
0135   * \implsparsesolverconcept
0136   *
0137   * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
0138   * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
0139   * and NumTraits<Scalar>::epsilon() for the tolerance.
0140   *
0141   * The tolerance corresponds to the relative residual error: |Ax-b|/|b|
0142   *
0143   * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format.
0144   * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled.
0145   * See \ref TopicMultiThreading for details.
0146   *
0147   * This class can be used as the direct solver classes. Here is a typical usage example:
0148   * \include BiCGSTAB_simple.cpp
0149   *
0150   * By default the iterations start with x=0 as an initial guess of the solution.
0151   * One can control the start using the solveWithGuess() method.
0152   *
0153   * BiCGSTAB can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
0154   *
0155   * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
0156   */
0157 template< typename _MatrixType, typename _Preconditioner>
0158 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
0159 {
0160   typedef IterativeSolverBase<BiCGSTAB> Base;
0161   using Base::matrix;
0162   using Base::m_error;
0163   using Base::m_iterations;
0164   using Base::m_info;
0165   using Base::m_isInitialized;
0166 public:
0167   typedef _MatrixType MatrixType;
0168   typedef typename MatrixType::Scalar Scalar;
0169   typedef typename MatrixType::RealScalar RealScalar;
0170   typedef _Preconditioner Preconditioner;
0171
0172 public:
0173
0174   /** Default constructor. */
0175   BiCGSTAB() : Base() {}
0176
0177   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
0178     *
0179     * This constructor is a shortcut for the default constructor followed
0180     * by a call to compute().
0181     *
0182     * \warning this class stores a reference to the matrix A as well as some
0183     * precomputed values that depend on it. Therefore, if \a A is changed
0184     * this class becomes invalid. Call compute() to update it with the new
0185     * matrix A, or modify a copy of A.
0186     */
0187   template<typename MatrixDerived>
0188   explicit BiCGSTAB(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
0189
0190   ~BiCGSTAB() {}
0191
0192   /** \internal */
0193   template<typename Rhs,typename Dest>
0194   void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
0195   {
0196     m_iterations = Base::maxIterations();
0197     m_error = Base::m_tolerance;
0198
0199     bool ret = internal::bicgstab(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error);
0200
0201     m_info = (!ret) ? NumericalIssue
0202            : m_error <= Base::m_tolerance ? Success
0203            : NoConvergence;
0204   }
0205
0206 protected:
0207
0208 };
0209
0210 } // end namespace RivetEigen
0211
0212 #endif // EIGEN_BICGSTAB_H