1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_BICGSTAB_H 12 #define EIGEN_BICGSTAB_H 13 14 namespace Eigen { 15 16 namespace internal { 17 18 /** \internal Low-level bi conjugate gradient stabilized algorithm 19 * \param mat The matrix A 20 * \param rhs The right hand side vector b 21 * \param x On input and initial solution, on output the computed solution. 22 * \param precond A preconditioner being able to efficiently solve for an 23 * approximation of Ax=b (regardless of b) 24 * \param iters On input the max number of iteration, on output the number of performed iterations. 25 * \param tol_error On input the tolerance error, on output an estimation of the relative error. 26 * \return false in the case of numerical issue, for example a break down of BiCGSTAB. 27 */ 28 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 29 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, 30 const Preconditioner& precond, Index& iters, 31 typename Dest::RealScalar& tol_error) 32 { 33 using std::sqrt; 34 using std::abs; 35 typedef typename Dest::RealScalar RealScalar; 36 typedef typename Dest::Scalar Scalar; 37 typedef Matrix<Scalar,Dynamic,1> VectorType; 38 RealScalar tol = tol_error; 39 Index maxIters = iters; 40 41 Index n = mat.cols(); 42 VectorType r = rhs - mat * x; 43 VectorType r0 = r; 44 45 RealScalar r0_sqnorm = r0.squaredNorm(); 46 RealScalar rhs_sqnorm = rhs.squaredNorm(); 47 if(rhs_sqnorm == 0) 48 { 49 x.setZero(); 50 return true; 51 } 52 Scalar rho = 1; 53 Scalar alpha = 1; 54 Scalar w = 1; 55 56 VectorType v = VectorType::Zero(n), p = VectorType::Zero(n); 57 VectorType y(n), z(n); 58 VectorType kt(n), ks(n); 59 60 VectorType s(n), t(n); 61 62 RealScalar tol2 = tol*tol*rhs_sqnorm; 63 RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); 64 Index i = 0; 65 Index restarts = 0; 66 67 while ( r.squaredNorm() > tol2 && i<maxIters ) 68 { 69 Scalar rho_old = rho; 70 71 rho = r0.dot(r); 72 if (abs(rho) < eps2*r0_sqnorm) 73 { 74 // The new residual vector became too orthogonal to the arbitrarily chosen direction r0 75 // Let's restart with a new r0: 76 r = rhs - mat * x; 77 r0 = r; 78 rho = r0_sqnorm = r.squaredNorm(); 79 if(restarts++ == 0) 80 i = 0; 81 } 82 Scalar beta = (rho/rho_old) * (alpha / w); 83 p = r + beta * (p - w * v); 84 85 y = precond.solve(p); 86 87 v.noalias() = mat * y; 88 89 alpha = rho / r0.dot(v); 90 s = r - alpha * v; 91 92 z = precond.solve(s); 93 t.noalias() = mat * z; 94 95 RealScalar tmp = t.squaredNorm(); 96 if(tmp>RealScalar(0)) 97 w = t.dot(s) / tmp; 98 else 99 w = Scalar(0); 100 x += alpha * y + w * z; 101 r = s - w * t; 102 ++i; 103 } 104 tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); 105 iters = i; 106 return true; 107 } 108 109 } 110 111 template< typename _MatrixType, 112 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 113 class BiCGSTAB; 114 115 namespace internal { 116 117 template< typename _MatrixType, typename _Preconditioner> 118 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> > 119 { 120 typedef _MatrixType MatrixType; 121 typedef _Preconditioner Preconditioner; 122 }; 123 124 } 125 126 /** \ingroup IterativeLinearSolvers_Module 127 * \brief A bi conjugate gradient stabilized solver for sparse square problems 128 * 129 * This class allows to solve for A.x = b sparse linear problems using a bi conjugate gradient 130 * stabilized algorithm. The vectors x and b can be either dense or sparse. 131 * 132 * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 133 * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 134 * 135 * \implsparsesolverconcept 136 * 137 * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 138 * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 139 * and NumTraits<Scalar>::epsilon() for the tolerance. 140 * 141 * The tolerance corresponds to the relative residual error: |Ax-b|/|b| 142 * 143 * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format. 144 * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled. 145 * See \ref TopicMultiThreading for details. 146 * 147 * This class can be used as the direct solver classes. Here is a typical usage example: 148 * \include BiCGSTAB_simple.cpp 149 * 150 * By default the iterations start with x=0 as an initial guess of the solution. 151 * One can control the start using the solveWithGuess() method. 152 * 153 * BiCGSTAB can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. 154 * 155 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 156 */ 157 template< typename _MatrixType, typename _Preconditioner> 158 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> > 159 { 160 typedef IterativeSolverBase<BiCGSTAB> Base; 161 using Base::matrix; 162 using Base::m_error; 163 using Base::m_iterations; 164 using Base::m_info; 165 using Base::m_isInitialized; 166 public: 167 typedef _MatrixType MatrixType; 168 typedef typename MatrixType::Scalar Scalar; 169 typedef typename MatrixType::RealScalar RealScalar; 170 typedef _Preconditioner Preconditioner; 171 172 public: 173 174 /** Default constructor. */ 175 BiCGSTAB() : Base() {} 176 177 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 178 * 179 * This constructor is a shortcut for the default constructor followed 180 * by a call to compute(). 181 * 182 * \warning this class stores a reference to the matrix A as well as some 183 * precomputed values that depend on it. Therefore, if \a A is changed 184 * this class becomes invalid. Call compute() to update it with the new 185 * matrix A, or modify a copy of A. 186 */ 187 template<typename MatrixDerived> 188 explicit BiCGSTAB(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 189 190 ~BiCGSTAB() {} 191 192 /** \internal */ 193 template<typename Rhs,typename Dest> 194 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 195 { 196 bool failed = false; 197 for(Index j=0; j<b.cols(); ++j) 198 { 199 m_iterations = Base::maxIterations(); 200 m_error = Base::m_tolerance; 201 202 typename Dest::ColXpr xj(x,j); 203 if(!internal::bicgstab(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) 204 failed = true; 205 } 206 m_info = failed ? NumericalIssue 207 : m_error <= Base::m_tolerance ? Success 208 : NoConvergence; 209 m_isInitialized = true; 210 } 211 212 /** \internal */ 213 using Base::_solve_impl; 214 template<typename Rhs,typename Dest> 215 void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const 216 { 217 x.resize(this->rows(),b.cols()); 218 x.setZero(); 219 _solve_with_guess_impl(b,x); 220 } 221 222 protected: 223 224 }; 225 226 } // end namespace Eigen 227 228 #endif // EIGEN_BICGSTAB_H 229
