$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2011 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, int& 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 int maxIters = iters; 00040 00041 int 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; 00063 RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); 00064 int i = 0; 00065 int restarts = 0; 00066 00067 while ( r.squaredNorm()/rhs_sqnorm > 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 choosen direction r0 00075 // Let's restart with a new r0: 00076 r0 = r; 00077 rho = r0_sqnorm = r.squaredNorm(); 00078 if(restarts++ == 0) 00079 i = 0; 00080 } 00081 Scalar beta = (rho/rho_old) * (alpha / w); 00082 p = r + beta * (p - w * v); 00083 00084 y = precond.solve(p); 00085 00086 v.noalias() = mat * y; 00087 00088 alpha = rho / r0.dot(v); 00089 s = r - alpha * v; 00090 00091 z = precond.solve(s); 00092 t.noalias() = mat * z; 00093 00094 RealScalar tmp = t.squaredNorm(); 00095 if(tmp>RealScalar(0)) 00096 w = t.dot(s) / tmp; 00097 else 00098 w = Scalar(0); 00099 x += alpha * y + w * z; 00100 r = s - w * t; 00101 ++i; 00102 } 00103 tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); 00104 iters = i; 00105 return true; 00106 } 00107 00108 } 00109 00110 template< typename _MatrixType, 00111 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 00112 class BiCGSTAB; 00113 00114 namespace internal { 00115 00116 template< typename _MatrixType, typename _Preconditioner> 00117 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> > 00118 { 00119 typedef _MatrixType MatrixType; 00120 typedef _Preconditioner Preconditioner; 00121 }; 00122 00123 } 00124 00158 template< typename _MatrixType, typename _Preconditioner> 00159 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> > 00160 { 00161 typedef IterativeSolverBase<BiCGSTAB> Base; 00162 using Base::mp_matrix; 00163 using Base::m_error; 00164 using Base::m_iterations; 00165 using Base::m_info; 00166 using Base::m_isInitialized; 00167 public: 00168 typedef _MatrixType MatrixType; 00169 typedef typename MatrixType::Scalar Scalar; 00170 typedef typename MatrixType::Index Index; 00171 typedef typename MatrixType::RealScalar RealScalar; 00172 typedef _Preconditioner Preconditioner; 00173 00174 public: 00175 00177 BiCGSTAB() : Base() {} 00178 00189 BiCGSTAB(const MatrixType& A) : Base(A) {} 00190 00191 ~BiCGSTAB() {} 00192 00198 template<typename Rhs,typename Guess> 00199 inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess> 00200 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const 00201 { 00202 eigen_assert(m_isInitialized && "BiCGSTAB is not initialized."); 00203 eigen_assert(Base::rows()==b.rows() 00204 && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b"); 00205 return internal::solve_retval_with_guess 00206 <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0); 00207 } 00208 00210 template<typename Rhs,typename Dest> 00211 void _solveWithGuess(const Rhs& b, Dest& x) const 00212 { 00213 bool failed = false; 00214 for(int j=0; j<b.cols(); ++j) 00215 { 00216 m_iterations = Base::maxIterations(); 00217 m_error = Base::m_tolerance; 00218 00219 typename Dest::ColXpr xj(x,j); 00220 if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) 00221 failed = true; 00222 } 00223 m_info = failed ? NumericalIssue 00224 : m_error <= Base::m_tolerance ? Success 00225 : NoConvergence; 00226 m_isInitialized = true; 00227 } 00228 00230 template<typename Rhs,typename Dest> 00231 void _solve(const Rhs& b, Dest& x) const 00232 { 00233 // x.setZero(); 00234 x = b; 00235 _solveWithGuess(b,x); 00236 } 00237 00238 protected: 00239 00240 }; 00241 00242 00243 namespace internal { 00244 00245 template<typename _MatrixType, typename _Preconditioner, typename Rhs> 00246 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> 00247 : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> 00248 { 00249 typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec; 00250 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00251 00252 template<typename Dest> void evalTo(Dest& dst) const 00253 { 00254 dec()._solve(rhs(),dst); 00255 } 00256 }; 00257 00258 } // end namespace internal 00259 00260 } // end namespace Eigen 00261 00262 #endif // EIGEN_BICGSTAB_H