BiCGSTAB.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
6 //
7 // Eigen is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Lesser General Public
9 // License as published by the Free Software Foundation; either
10 // version 3 of the License, or (at your option) any later version.
11 //
12 // Alternatively, you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License as
14 // published by the Free Software Foundation; either version 2 of
15 // the License, or (at your option) any later version.
16 //
17 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License and a copy of the GNU General Public License along with
24 // Eigen. If not, see <http://www.gnu.org/licenses/>.
25 
26 #ifndef EIGEN_BICGSTAB_H
27 #define EIGEN_BICGSTAB_H
28 
29 namespace Eigen {
30 
31 namespace internal {
32 
43 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
44 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
45  const Preconditioner& precond, int& iters,
46  typename Dest::RealScalar& tol_error)
47 {
48  using std::sqrt;
49  using std::abs;
50  typedef typename Dest::RealScalar RealScalar;
51  typedef typename Dest::Scalar Scalar;
52  typedef Matrix<Scalar,Dynamic,1> VectorType;
53  RealScalar tol = tol_error;
54  int maxIters = iters;
55 
56  int n = mat.cols();
57  VectorType r = rhs - mat * x;
58  VectorType r0 = r;
59 
60  RealScalar r0_sqnorm = r0.squaredNorm();
61  Scalar rho = 1;
62  Scalar alpha = 1;
63  Scalar w = 1;
64 
65  VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
66  VectorType y(n), z(n);
67  VectorType kt(n), ks(n);
68 
69  VectorType s(n), t(n);
70 
71  RealScalar tol2 = tol*tol;
72  int i = 0;
73 
74  while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
75  {
76  Scalar rho_old = rho;
77 
78  rho = r0.dot(r);
79  if (rho == Scalar(0)) return false; /* New search directions cannot be found */
80  Scalar beta = (rho/rho_old) * (alpha / w);
81  p = r + beta * (p - w * v);
82 
83  y = precond.solve(p);
84 
85  v.noalias() = mat * y;
86 
87  alpha = rho / r0.dot(v);
88  s = r - alpha * v;
89 
90  z = precond.solve(s);
91  t.noalias() = mat * z;
92 
93  w = t.dot(s) / t.squaredNorm();
94  x += alpha * y + w * z;
95  r = s - w * t;
96  ++i;
97  }
98  tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
99  iters = i;
100  return true;
101 }
102 
103 }
104 
105 template< typename _MatrixType,
106  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
107 class BiCGSTAB;
108 
109 namespace internal {
110 
111 template< typename _MatrixType, typename _Preconditioner>
112 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
113 {
114  typedef _MatrixType MatrixType;
115  typedef _Preconditioner Preconditioner;
116 };
117 
118 }
119 
166 template< typename _MatrixType, typename _Preconditioner>
167 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
168 {
170  using Base::mp_matrix;
171  using Base::m_error;
172  using Base::m_iterations;
173  using Base::m_info;
174  using Base::m_isInitialized;
175 public:
176  typedef _MatrixType MatrixType;
177  typedef typename MatrixType::Scalar Scalar;
178  typedef typename MatrixType::Index Index;
179  typedef typename MatrixType::RealScalar RealScalar;
180  typedef _Preconditioner Preconditioner;
181 
182 public:
183 
185  BiCGSTAB() : Base() {}
186 
197  BiCGSTAB(const MatrixType& A) : Base(A) {}
198 
200 
206  template<typename Rhs,typename Guess>
207  inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
208  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
209  {
210  eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
211  eigen_assert(Base::rows()==b.rows()
212  && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
213  return internal::solve_retval_with_guess
214  <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
215  }
216 
218  template<typename Rhs,typename Dest>
219  void _solveWithGuess(const Rhs& b, Dest& x) const
220  {
221  bool failed = false;
222  for(int j=0; j<b.cols(); ++j)
223  {
226 
227  typename Dest::ColXpr xj(x,j);
229  failed = true;
230  }
231  m_info = failed ? NumericalIssue
233  : NoConvergence;
234  m_isInitialized = true;
235  }
236 
238  template<typename Rhs,typename Dest>
239  void _solve(const Rhs& b, Dest& x) const
240  {
241  x.setZero();
242  _solveWithGuess(b,x);
243  }
244 
245 protected:
246 
247 };
248 
249 
250 namespace internal {
251 
252  template<typename _MatrixType, typename _Preconditioner, typename Rhs>
253 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
254  : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
255 {
256  typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
257  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
258 
259  template<typename Dest> void evalTo(Dest& dst) const
260  {
261  dec()._solve(rhs(),dst);
262  }
263 };
264 
265 } // end namespace internal
266 
267 } // end namespace Eigen
268 
269 #endif // EIGEN_BICGSTAB_H