Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
BiCGSTAB.h
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@inria.fr>
5 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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 #include "./InternalHeaderCheck.h"
15 
16 namespace Eigen {
17 
18 namespace internal {
19 
30 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
31 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
32  const Preconditioner& precond, Index& iters,
33  typename Dest::RealScalar& tol_error)
34 {
35  using std::sqrt;
36  using std::abs;
37  typedef typename Dest::RealScalar RealScalar;
38  typedef typename Dest::Scalar Scalar;
39  typedef Matrix<Scalar,Dynamic,1> VectorType;
40  RealScalar tol = tol_error;
41  Index maxIters = iters;
42 
43  Index n = mat.cols();
44  VectorType r = rhs - mat * x;
45  VectorType r0 = r;
46 
47  RealScalar r0_sqnorm = r0.squaredNorm();
48  RealScalar rhs_sqnorm = rhs.squaredNorm();
49  if(rhs_sqnorm == 0)
50  {
51  x.setZero();
52  return true;
53  }
54  Scalar rho = 1;
55  Scalar alpha = 1;
56  Scalar w = 1;
57 
58  VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
59  VectorType y(n), z(n);
60  VectorType kt(n), ks(n);
61 
62  VectorType s(n), t(n);
63 
64  RealScalar tol2 = tol*tol*rhs_sqnorm;
65  RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
66  Index i = 0;
67  Index restarts = 0;
68 
69  while ( r.squaredNorm() > tol2 && i<maxIters )
70  {
71  Scalar rho_old = rho;
72 
73  rho = r0.dot(r);
74  if (abs(rho) < eps2*r0_sqnorm)
75  {
76  // The new residual vector became too orthogonal to the arbitrarily chosen direction r0
77  // Let's restart with a new r0:
78  r = rhs - mat * x;
79  r0 = r;
80  rho = r0_sqnorm = r.squaredNorm();
81  if(restarts++ == 0)
82  i = 0;
83  }
84  Scalar beta = (rho/rho_old) * (alpha / w);
85  p = r + beta * (p - w * v);
86 
87  y = precond.solve(p);
88 
89  v.noalias() = mat * y;
90 
91  alpha = rho / r0.dot(v);
92  s = r - alpha * v;
93 
94  z = precond.solve(s);
95  t.noalias() = mat * z;
96 
97  RealScalar tmp = t.squaredNorm();
98  if(tmp>RealScalar(0))
99  w = t.dot(s) / tmp;
100  else
101  w = Scalar(0);
102  x += alpha * y + w * z;
103  r = s - w * t;
104  ++i;
105  }
106  tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
107  iters = i;
108  return true;
109 }
110 
111 }
112 
113 template< typename MatrixType_,
114  typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
115 class BiCGSTAB;
116 
117 namespace internal {
118 
119 template< typename MatrixType_, typename Preconditioner_>
120 struct traits<BiCGSTAB<MatrixType_,Preconditioner_> >
121 {
122  typedef MatrixType_ MatrixType;
123  typedef Preconditioner_ Preconditioner;
124 };
125 
126 }
127 
159 template< typename MatrixType_, typename Preconditioner_>
160 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<MatrixType_,Preconditioner_> >
161 {
163  using Base::matrix;
164  using Base::m_error;
165  using Base::m_iterations;
166  using Base::m_info;
167  using Base::m_isInitialized;
168 public:
169  typedef MatrixType_ MatrixType;
170  typedef typename MatrixType::Scalar Scalar;
171  typedef typename MatrixType::RealScalar RealScalar;
172  typedef Preconditioner_ Preconditioner;
173 
174 public:
175 
177  BiCGSTAB() : Base() {}
178 
189  template<typename MatrixDerived>
190  explicit BiCGSTAB(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
191 
192  ~BiCGSTAB() {}
193 
195  template<typename Rhs,typename Dest>
196  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
197  {
198  m_iterations = Base::maxIterations();
199  m_error = Base::m_tolerance;
200 
201  bool ret = internal::bicgstab(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error);
202 
203  m_info = (!ret) ? NumericalIssue
204  : m_error <= Base::m_tolerance ? Success
205  : NoConvergence;
206  }
207 
208 protected:
209 
210 };
211 
212 } // end namespace Eigen
213 
214 #endif // EIGEN_BICGSTAB_H
A bi conjugate gradient stabilized solver for sparse square problems.
Definition: BiCGSTAB.h:161
BiCGSTAB(const EigenBase< MatrixDerived > &A)
Definition: BiCGSTAB.h:190
BiCGSTAB()
Definition: BiCGSTAB.h:177
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:146
Index maxIterations() const
Definition: IterativeSolverBase.h:286
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
@ NoConvergence
Definition: Constants.h:448
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32