![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_BICGSTAB_H 00012 #define EIGEN_BICGSTAB_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00028 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 00029 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, 00030 const Preconditioner& precond, Index& iters, 00031 typename Dest::RealScalar& tol_error) 00032 { 00033 using std::sqrt; 00034 using std::abs; 00035 typedef typename Dest::RealScalar RealScalar; 00036 typedef typename Dest::Scalar Scalar; 00037 typedef Matrix<Scalar,Dynamic,1> VectorType; 00038 RealScalar tol = tol_error; 00039 Index maxIters = iters; 00040 00041 Index n = mat.cols(); 00042 VectorType r = rhs - mat * x; 00043 VectorType r0 = r; 00044 00045 RealScalar r0_sqnorm = r0.squaredNorm(); 00046 RealScalar rhs_sqnorm = rhs.squaredNorm(); 00047 if(rhs_sqnorm == 0) 00048 { 00049 x.setZero(); 00050 return true; 00051 } 00052 Scalar rho = 1; 00053 Scalar alpha = 1; 00054 Scalar w = 1; 00055 00056 VectorType v = VectorType::Zero(n), p = VectorType::Zero(n); 00057 VectorType y(n), z(n); 00058 VectorType kt(n), ks(n); 00059 00060 VectorType s(n), t(n); 00061 00062 RealScalar tol2 = tol*tol*rhs_sqnorm; 00063 RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); 00064 Index i = 0; 00065 Index restarts = 0; 00066 00067 while ( r.squaredNorm() > tol2 && i<maxIters ) 00068 { 00069 Scalar rho_old = rho; 00070 00071 rho = r0.dot(r); 00072 if (abs(rho) < eps2*r0_sqnorm) 00073 { 00074 // The new residual vector became too orthogonal to the arbitrarily chosen direction r0 00075 // Let's restart with a new r0: 00076 r = rhs - mat * x; 00077 r0 = r; 00078 rho = r0_sqnorm = r.squaredNorm(); 00079 if(restarts++ == 0) 00080 i = 0; 00081 } 00082 Scalar beta = (rho/rho_old) * (alpha / w); 00083 p = r + beta * (p - w * v); 00084 00085 y = precond.solve(p); 00086 00087 v.noalias() = mat * y; 00088 00089 alpha = rho / r0.dot(v); 00090 s = r - alpha * v; 00091 00092 z = precond.solve(s); 00093 t.noalias() = mat * z; 00094 00095 RealScalar tmp = t.squaredNorm(); 00096 if(tmp>RealScalar(0)) 00097 w = t.dot(s) / tmp; 00098 else 00099 w = Scalar(0); 00100 x += alpha * y + w * z; 00101 r = s - w * t; 00102 ++i; 00103 } 00104 tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); 00105 iters = i; 00106 return true; 00107 } 00108 00109 } 00110 00111 template< typename _MatrixType, 00112 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 00113 class BiCGSTAB; 00114 00115 namespace internal { 00116 00117 template< typename _MatrixType, typename _Preconditioner> 00118 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> > 00119 { 00120 typedef _MatrixType MatrixType; 00121 typedef _Preconditioner Preconditioner; 00122 }; 00123 00124 } 00125 00157 template< typename _MatrixType, typename _Preconditioner> 00158 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> > 00159 { 00160 typedef IterativeSolverBase<BiCGSTAB> Base; 00161 using Base::matrix; 00162 using Base::m_error; 00163 using Base::m_iterations; 00164 using Base::m_info; 00165 using Base::m_isInitialized; 00166 public: 00167 typedef _MatrixType MatrixType; 00168 typedef typename MatrixType::Scalar Scalar; 00169 typedef typename MatrixType::RealScalar RealScalar; 00170 typedef _Preconditioner Preconditioner; 00171 00172 public: 00173 00175 BiCGSTAB() : Base() {} 00176 00187 template<typename MatrixDerived> 00188 explicit BiCGSTAB(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 00189 00190 ~BiCGSTAB() {} 00191 00193 template<typename Rhs,typename Dest> 00194 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 00195 { 00196 bool failed = false; 00197 for(Index j=0; j<b.cols(); ++j) 00198 { 00199 m_iterations = Base::maxIterations(); 00200 m_error = Base::m_tolerance; 00201 00202 typename Dest::ColXpr xj(x,j); 00203 if(!internal::bicgstab(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) 00204 failed = true; 00205 } 00206 m_info = failed ? NumericalIssue 00207 : m_error <= Base::m_tolerance ? Success 00208 : NoConvergence; 00209 m_isInitialized = true; 00210 } 00211 00213 using Base::_solve_impl; 00214 template<typename Rhs,typename Dest> 00215 void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const 00216 { 00217 x.resize(this->rows(),b.cols()); 00218 x.setZero(); 00219 _solve_with_guess_impl(b,x); 00220 } 00221 00222 protected: 00223 00224 }; 00225 00226 } // end namespace Eigen 00227 00228 #endif // EIGEN_BICGSTAB_H