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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_BICGSTAB_H
00027 #define EIGEN_BICGSTAB_H
00028 
00029 namespace Eigen { 
00030 
00031 namespace internal {
00032 
00043 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
00044 bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
00045               const Preconditioner& precond, int& iters,
00046               typename Dest::RealScalar& tol_error)
00047 {
00048   using std::sqrt;
00049   using std::abs;
00050   typedef typename Dest::RealScalar RealScalar;
00051   typedef typename Dest::Scalar Scalar;
00052   typedef Matrix<Scalar,Dynamic,1> VectorType;
00053   RealScalar tol = tol_error;
00054   int maxIters = iters;
00055 
00056   int n = mat.cols();
00057   VectorType r  = rhs - mat * x;
00058   VectorType r0 = r;
00059   
00060   RealScalar r0_sqnorm = r0.squaredNorm();
00061   Scalar rho    = 1;
00062   Scalar alpha  = 1;
00063   Scalar w      = 1;
00064   
00065   VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
00066   VectorType y(n),  z(n);
00067   VectorType kt(n), ks(n);
00068 
00069   VectorType s(n), t(n);
00070 
00071   RealScalar tol2 = tol*tol;
00072   int i = 0;
00073 
00074   while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
00075   {
00076     Scalar rho_old = rho;
00077 
00078     rho = r0.dot(r);
00079     if (rho == Scalar(0)) return false; /* New search directions cannot be found */
00080     Scalar beta = (rho/rho_old) * (alpha / w);
00081     p = r + beta * (p - w * v);
00082     
00083     y = precond.solve(p);
00084     
00085     v.noalias() = mat * y;
00086 
00087     alpha = rho / r0.dot(v);
00088     s = r - alpha * v;
00089 
00090     z = precond.solve(s);
00091     t.noalias() = mat * z;
00092 
00093     w = t.dot(s) / t.squaredNorm();
00094     x += alpha * y + w * z;
00095     r = s - w * t;
00096     ++i;
00097   }
00098   tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
00099   iters = i;
00100   return true; 
00101 }
00102 
00103 }
00104 
00105 template< typename _MatrixType,
00106           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
00107 class BiCGSTAB;
00108 
00109 namespace internal {
00110 
00111 template< typename _MatrixType, typename _Preconditioner>
00112 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
00113 {
00114   typedef _MatrixType MatrixType;
00115   typedef _Preconditioner Preconditioner;
00116 };
00117 
00118 }
00119 
00166 template< typename _MatrixType, typename _Preconditioner>
00167 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
00168 {
00169   typedef IterativeSolverBase<BiCGSTAB> Base;
00170   using Base::mp_matrix;
00171   using Base::m_error;
00172   using Base::m_iterations;
00173   using Base::m_info;
00174   using Base::m_isInitialized;
00175 public:
00176   typedef _MatrixType MatrixType;
00177   typedef typename MatrixType::Scalar Scalar;
00178   typedef typename MatrixType::Index Index;
00179   typedef typename MatrixType::RealScalar RealScalar;
00180   typedef _Preconditioner Preconditioner;
00181 
00182 public:
00183 
00185   BiCGSTAB() : Base() {}
00186 
00197   BiCGSTAB(const MatrixType& A) : Base(A) {}
00198 
00199   ~BiCGSTAB() {}
00200   
00206   template<typename Rhs,typename Guess>
00207   inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
00208   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00209   {
00210     eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
00211     eigen_assert(Base::rows()==b.rows()
00212               && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
00213     return internal::solve_retval_with_guess
00214             <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
00215   }
00216   
00218   template<typename Rhs,typename Dest>
00219   void _solveWithGuess(const Rhs& b, Dest& x) const
00220   {    
00221     bool failed = false;
00222     for(int j=0; j<b.cols(); ++j)
00223     {
00224       m_iterations = Base::maxIterations();
00225       m_error = Base::m_tolerance;
00226       
00227       typename Dest::ColXpr xj(x,j);
00228       if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error))
00229         failed = true;
00230     }
00231     m_info = failed ? NumericalIssue
00232            : m_error <= Base::m_tolerance ? Success
00233            : NoConvergence;
00234     m_isInitialized = true;
00235   }
00236 
00238   template<typename Rhs,typename Dest>
00239   void _solve(const Rhs& b, Dest& x) const
00240   {
00241     x.setZero();
00242     _solveWithGuess(b,x);
00243   }
00244 
00245 protected:
00246 
00247 };
00248 
00249 
00250 namespace internal {
00251 
00252   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
00253 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00254   : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00255 {
00256   typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
00257   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00258 
00259   template<typename Dest> void evalTo(Dest& dst) const
00260   {
00261     dec()._solve(rhs(),dst);
00262   }
00263 };
00264 
00265 } // end namespace internal
00266 
00267 } // end namespace Eigen
00268 
00269 #endif // EIGEN_BICGSTAB_H