00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
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;
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 }
00266
00267 }
00268
00269 #endif // EIGEN_BICGSTAB_H