11 #ifndef EIGEN_BICGSTAB_H
12 #define EIGEN_BICGSTAB_H
14 #include "./InternalHeaderCheck.h"
30 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
31 bool bicgstab(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
32 const Preconditioner& precond,
Index& iters,
33 typename Dest::RealScalar& tol_error)
37 typedef typename Dest::RealScalar RealScalar;
38 typedef typename Dest::Scalar Scalar;
39 typedef Matrix<Scalar,Dynamic,1> VectorType;
40 RealScalar tol = tol_error;
41 Index maxIters = iters;
44 VectorType r = rhs - mat * x;
47 RealScalar r0_sqnorm = r0.squaredNorm();
48 RealScalar rhs_sqnorm = rhs.squaredNorm();
58 VectorType v = VectorType::Zero(n), p = VectorType::Zero(n);
59 VectorType y(n), z(n);
60 VectorType kt(n), ks(n);
62 VectorType s(n), t(n);
64 RealScalar tol2 = tol*tol*rhs_sqnorm;
65 RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
69 while ( r.squaredNorm() > tol2 && i<maxIters )
74 if (
abs(rho) < eps2*r0_sqnorm)
80 rho = r0_sqnorm = r.squaredNorm();
84 Scalar beta = (rho/rho_old) * (alpha / w);
85 p = r + beta * (p - w * v);
89 v.noalias() = mat * y;
91 alpha = rho / r0.dot(v);
95 t.noalias() = mat * z;
97 RealScalar tmp = t.squaredNorm();
102 x += alpha * y + w * z;
106 tol_error =
sqrt(r.squaredNorm()/rhs_sqnorm);
113 template<
typename MatrixType_,
114 typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
119 template<
typename MatrixType_,
typename Preconditioner_>
120 struct traits<BiCGSTAB<MatrixType_,Preconditioner_> >
122 typedef MatrixType_ MatrixType;
123 typedef Preconditioner_ Preconditioner;
159 template<
typename MatrixType_,
typename Preconditioner_>
165 using Base::m_iterations;
167 using Base::m_isInitialized;
169 typedef MatrixType_ MatrixType;
170 typedef typename MatrixType::Scalar Scalar;
171 typedef typename MatrixType::RealScalar RealScalar;
172 typedef Preconditioner_ Preconditioner;
189 template<
typename MatrixDerived>
195 template<
typename Rhs,
typename Dest>
196 void _solve_vector_with_guess_impl(
const Rhs& b, Dest& x)
const
199 m_error = Base::m_tolerance;
201 bool ret = internal::bicgstab(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error);
204 : m_error <= Base::m_tolerance ?
Success
A bi conjugate gradient stabilized solver for sparse square problems.
Definition: BiCGSTAB.h:161
BiCGSTAB(const EigenBase< MatrixDerived > &A)
Definition: BiCGSTAB.h:190
BiCGSTAB()
Definition: BiCGSTAB.h:177
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:146
Index maxIterations() const
Definition: IterativeSolverBase.h:286
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
@ NoConvergence
Definition: Constants.h:448
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32