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   VectorType r  = rhs - mat * x;
00043   VectorType r0 = r;
00044   
00045   RealScalar r0_sqnorm = r0.squaredNorm();
00046   Scalar rho    = 1;
00047   Scalar alpha  = 1;
00048   Scalar w      = 1;
00049   
00050   VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
00051   VectorType y(n),  z(n);
00052   VectorType kt(n), ks(n);
00053 
00054   VectorType s(n), t(n);
00055 
00056   RealScalar tol2 = tol*tol;
00057   int i = 0;
00058 
00059   while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
00060   {
00061     Scalar rho_old = rho;
00062 
00063     rho = r0.dot(r);
00064     if (rho == Scalar(0)) return false; /* New search directions cannot be found */
00065     Scalar beta = (rho/rho_old) * (alpha / w);
00066     p = r + beta * (p - w * v);
00067     
00068     y = precond.solve(p);
00069     
00070     v.noalias() = mat * y;
00071 
00072     alpha = rho / r0.dot(v);
00073     s = r - alpha * v;
00074 
00075     z = precond.solve(s);
00076     t.noalias() = mat * z;
00077 
00078     w = t.dot(s) / t.squaredNorm();
00079     x += alpha * y + w * z;
00080     r = s - w * t;
00081     ++i;
00082   }
00083   tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
00084   iters = i;
00085   return true; 
00086 }
00087 
00088 }
00089 
00090 template< typename _MatrixType,
00091           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
00092 class BiCGSTAB;
00093 
00094 namespace internal {
00095 
00096 template< typename _MatrixType, typename _Preconditioner>
00097 struct traits<BiCGSTAB<_MatrixType,_Preconditioner> >
00098 {
00099   typedef _MatrixType MatrixType;
00100   typedef _Preconditioner Preconditioner;
00101 };
00102 
00103 }
00104 
00151 template< typename _MatrixType, typename _Preconditioner>
00152 class BiCGSTAB : public IterativeSolverBase<BiCGSTAB<_MatrixType,_Preconditioner> >
00153 {
00154   typedef IterativeSolverBase<BiCGSTAB> Base;
00155   using Base::mp_matrix;
00156   using Base::m_error;
00157   using Base::m_iterations;
00158   using Base::m_info;
00159   using Base::m_isInitialized;
00160 public:
00161   typedef _MatrixType MatrixType;
00162   typedef typename MatrixType::Scalar Scalar;
00163   typedef typename MatrixType::Index Index;
00164   typedef typename MatrixType::RealScalar RealScalar;
00165   typedef _Preconditioner Preconditioner;
00166 
00167 public:
00168 
00170   BiCGSTAB() : Base() {}
00171 
00182   BiCGSTAB(const MatrixType& A) : Base(A) {}
00183 
00184   ~BiCGSTAB() {}
00185   
00191   template<typename Rhs,typename Guess>
00192   inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess>
00193   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00194   {
00195     eigen_assert(m_isInitialized && "BiCGSTAB is not initialized.");
00196     eigen_assert(Base::rows()==b.rows()
00197               && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b");
00198     return internal::solve_retval_with_guess
00199             <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0);
00200   }
00201   
00203   template<typename Rhs,typename Dest>
00204   void _solveWithGuess(const Rhs& b, Dest& x) const
00205   {    
00206     bool failed = false;
00207     for(int j=0; j<b.cols(); ++j)
00208     {
00209       m_iterations = Base::maxIterations();
00210       m_error = Base::m_tolerance;
00211       
00212       typename Dest::ColXpr xj(x,j);
00213       if(!internal::bicgstab(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error))
00214         failed = true;
00215     }
00216     m_info = failed ? NumericalIssue
00217            : m_error <= Base::m_tolerance ? Success
00218            : NoConvergence;
00219     m_isInitialized = true;
00220   }
00221 
00223   template<typename Rhs,typename Dest>
00224   void _solve(const Rhs& b, Dest& x) const
00225   {
00226     x.setZero();
00227     _solveWithGuess(b,x);
00228   }
00229 
00230 protected:
00231 
00232 };
00233 
00234 
00235 namespace internal {
00236 
00237   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
00238 struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00239   : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs>
00240 {
00241   typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec;
00242   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00243 
00244   template<typename Dest> void evalTo(Dest& dst) const
00245   {
00246     dec()._solve(rhs(),dst);
00247   }
00248 };
00249 
00250 } // end namespace internal
00251 
00252 } // end namespace Eigen
00253 
00254 #endif // EIGEN_BICGSTAB_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:12