BiCGSTAB.h
Go to the documentation of this file.
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   x = precond.solve(x);
00043   VectorType r  = rhs - mat * x;
00044   VectorType r0 = r;
00045   
00046   RealScalar r0_sqnorm = r0.squaredNorm();
00047   RealScalar rhs_sqnorm = rhs.squaredNorm();
00048   if(rhs_sqnorm == 0)
00049   {
00050     x.setZero();
00051     return true;
00052   }
00053   Scalar rho    = 1;
00054   Scalar alpha  = 1;
00055   Scalar w      = 1;
00056   
00057   VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
00058   VectorType y(n),  z(n);
00059   VectorType kt(n), ks(n);
00060 
00061   VectorType s(n), t(n);
00062 
00063   RealScalar tol2 = tol*tol;
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 (internal::isMuchSmallerThan(rho,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 
00171 template< typename _MatrixType, typename _Preconditioner>
00172 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
00173 {
00174   typedef IterativeSolverBase<BiCGSTAB> Base;
00175   using Base::mp_matrix;
00176   using Base::m_error;
00177   using Base::m_iterations;
00178   using Base::m_info;
00179   using Base::m_isInitialized;
00180 public:
00181   typedef _MatrixType MatrixType;
00182   typedef typename MatrixType::Scalar Scalar;
00183   typedef typename MatrixType::Index Index;
00184   typedef typename MatrixType::RealScalar RealScalar;
00185   typedef _Preconditioner Preconditioner;
00186 
00187 public:
00188 
00190   BiCGSTAB() : Base() {}
00191 
00202   BiCGSTAB(const MatrixType& A) : Base(A) {}
00203 
00204   ~BiCGSTAB() {}
00205   
00211   template<typename Rhs,typename Guess>
00212   inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
00213   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00214   {
00215     eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
00216     eigen_assert(Base::rows()==b.rows()
00217               && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
00218     return internal::solve_retval_with_guess
00219             <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
00220   }
00221   
00223   template<typename Rhs,typename Dest>
00224   void _solveWithGuess(const Rhs& b, Dest& x) const
00225   {    
00226     bool failed = false;
00227     for(int j=0; j<b.cols(); ++j)
00228     {
00229       m_iterations = Base::maxIterations();
00230       m_error = Base::m_tolerance;
00231       
00232       typename Dest::ColXpr xj(x,j);
00233       if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error))
00234         failed = true;
00235     }
00236     m_info = failed ? NumericalIssue
00237            : m_error <= Base::m_tolerance ? Success
00238            : NoConvergence;
00239     m_isInitialized = true;
00240   }
00241 
00243   template<typename Rhs,typename Dest>
00244   void _solve(const Rhs& b, Dest& x) const
00245   {
00246 //     x.setZero();
00247   x = b;
00248     _solveWithGuess(b,x);
00249   }
00250 
00251 protected:
00252 
00253 };
00254 
00255 
00256 namespace internal {
00257 
00258   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
00259 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00260   : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00261 {
00262   typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
00263   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00264 
00265   template<typename Dest> void evalTo(Dest& dst) const
00266   {
00267     dec()._solve(rhs(),dst);
00268   }
00269 };
00270 
00271 } // end namespace internal
00272 
00273 } // end namespace Eigen
00274 
00275 #endif // EIGEN_BICGSTAB_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:57:53