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 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,Index & iters,typename Dest::RealScalar & tol_error)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