Eigen  3.3.3
BiCGSTAB.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends