Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
16namespace Eigen {
17
18namespace internal {
19
30template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
31bool 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
113template< typename MatrixType_,
114 typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
115class BiCGSTAB;
116
117namespace internal {
118
119template< typename MatrixType_, typename Preconditioner_>
120struct traits<BiCGSTAB<MatrixType_,Preconditioner_> >
121{
122 typedef MatrixType_ MatrixType;
123 typedef Preconditioner_ Preconditioner;
124};
125
126}
127
159template< typename MatrixType_, typename Preconditioner_>
160class 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;
168public:
169 typedef MatrixType_ MatrixType;
170 typedef typename MatrixType::Scalar Scalar;
171 typedef typename MatrixType::RealScalar RealScalar;
172 typedef Preconditioner_ Preconditioner;
173
174public:
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
206 }
207
208protected:
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:283
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
@ NoConvergence
Definition: Constants.h:448
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32