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 // 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>
bicgstab(const MatrixType & mat,const Rhs & rhs,Dest & x,const Preconditioner & precond,int & iters,typename Dest::RealScalar & tol_error)29 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
30               const Preconditioner& precond, int& 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   int maxIters = iters;
40 
41   int 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;
63   RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
64   int i = 0;
65   int restarts = 0;
66 
67   while ( r.squaredNorm()/rhs_sqnorm > 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 choosen direction r0
75       // Let's restart with a new r0:
76       r0 = r;
77       rho = r0_sqnorm = r.squaredNorm();
78       if(restarts++ == 0)
79         i = 0;
80     }
81     Scalar beta = (rho/rho_old) * (alpha / w);
82     p = r + beta * (p - w * v);
83 
84     y = precond.solve(p);
85 
86     v.noalias() = mat * y;
87 
88     alpha = rho / r0.dot(v);
89     s = r - alpha * v;
90 
91     z = precond.solve(s);
92     t.noalias() = mat * z;
93 
94     RealScalar tmp = t.squaredNorm();
95     if(tmp>RealScalar(0))
96       w = t.dot(s) / tmp;
97     else
98       w = Scalar(0);
99     x += alpha * y + w * z;
100     r = s - w * t;
101     ++i;
102   }
103   tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
104   iters = i;
105   return true;
106 }
107 
108 }
109 
110 template< typename _MatrixType,
111           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
112 class BiCGSTAB;
113 
114 namespace internal {
115 
116 template< typename _MatrixType, typename _Preconditioner>
117 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
118 {
119   typedef _MatrixType MatrixType;
120   typedef _Preconditioner Preconditioner;
121 };
122 
123 }
124 
125 /** \ingroup IterativeLinearSolvers_Module
126   * \brief A bi conjugate gradient stabilized solver for sparse square problems
127   *
128   * This class allows to solve for A.x = b sparse linear problems using a bi conjugate gradient
129   * stabilized algorithm. The vectors x and b can be either dense or sparse.
130   *
131   * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
132   * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
133   *
134   * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
135   * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
136   * and NumTraits<Scalar>::epsilon() for the tolerance.
137   *
138   * This class can be used as the direct solver classes. Here is a typical usage example:
139   * \code
140   * int n = 10000;
141   * VectorXd x(n), b(n);
142   * SparseMatrix<double> A(n,n);
143   * // fill A and b
144   * BiCGSTAB<SparseMatrix<double> > solver;
145   * solver.compute(A);
146   * x = solver.solve(b);
147   * std::cout << "#iterations:     " << solver.iterations() << std::endl;
148   * std::cout << "estimated error: " << solver.error()      << std::endl;
149   * // update b, and solve again
150   * x = solver.solve(b);
151   * \endcode
152   *
153   * By default the iterations start with x=0 as an initial guess of the solution.
154   * One can control the start using the solveWithGuess() method. Here is a step by
155   * step execution example starting with a random guess and printing the evolution
156   * of the estimated error:
157   * * \code
158   * x = VectorXd::Random(n);
159   * solver.setMaxIterations(1);
160   * int i = 0;
161   * do {
162   *   x = solver.solveWithGuess(b,x);
163   *   std::cout << i << " : " << solver.error() << std::endl;
164   *   ++i;
165   * } while (solver.info()!=Success && i<100);
166   * \endcode
167   * Note that such a step by step excution is slightly slower.
168   *
169   * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
170   */
171 template< typename _MatrixType, typename _Preconditioner>
172 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
173 {
174   typedef IterativeSolverBase<BiCGSTAB> Base;
175   using Base::mp_matrix;
176   using Base::m_error;
177   using Base::m_iterations;
178   using Base::m_info;
179   using Base::m_isInitialized;
180 public:
181   typedef _MatrixType MatrixType;
182   typedef typename MatrixType::Scalar Scalar;
183   typedef typename MatrixType::Index Index;
184   typedef typename MatrixType::RealScalar RealScalar;
185   typedef _Preconditioner Preconditioner;
186 
187 public:
188 
189   /** Default constructor. */
190   BiCGSTAB() : Base() {}
191 
192   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
193     *
194     * This constructor is a shortcut for the default constructor followed
195     * by a call to compute().
196     *
197     * \warning this class stores a reference to the matrix A as well as some
198     * precomputed values that depend on it. Therefore, if \a A is changed
199     * this class becomes invalid. Call compute() to update it with the new
200     * matrix A, or modify a copy of A.
201     */
202   BiCGSTAB(const MatrixType& A) : Base(A) {}
203 
204   ~BiCGSTAB() {}
205 
206   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
207     * \a x0 as an initial solution.
208     *
209     * \sa compute()
210     */
211   template<typename Rhs,typename Guess>
212   inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
213   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
214   {
215     eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
216     eigen_assert(Base::rows()==b.rows()
217               && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
218     return internal::solve_retval_with_guess
219             <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
220   }
221 
222   /** \internal */
223   template<typename Rhs,typename Dest>
224   void _solveWithGuess(const Rhs& b, Dest& x) const
225   {
226     bool failed = false;
227     for(int j=0; j<b.cols(); ++j)
228     {
229       m_iterations = Base::maxIterations();
230       m_error = Base::m_tolerance;
231 
232       typename Dest::ColXpr xj(x,j);
233       if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error))
234         failed = true;
235     }
236     m_info = failed ? NumericalIssue
237            : m_error <= Base::m_tolerance ? Success
238            : NoConvergence;
239     m_isInitialized = true;
240   }
241 
242   /** \internal */
243   template<typename Rhs,typename Dest>
244   void _solve(const Rhs& b, Dest& x) const
245   {
246 //     x.setZero();
247   x = b;
248     _solveWithGuess(b,x);
249   }
250 
251 protected:
252 
253 };
254 
255 
256 namespace internal {
257 
258   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
259 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
260   : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
261 {
262   typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
263   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
264 
265   template<typename Dest> void evalTo(Dest& dst) const
266   {
267     dec()._solve(rhs(),dst);
268   }
269 };
270 
271 } // end namespace internal
272 
273 } // end namespace Eigen
274 
275 #endif // EIGEN_BICGSTAB_H
276