1 #ifndef VIENNACL_LINALG_BICGSTAB_HPP_
2 #define VIENNACL_LINALG_BICGSTAB_HPP_
53 : tol_(tol), iterations_(max_iters), iterations_before_restart_(max_iters_before_restart) {}
67 double error()
const {
return last_error_; }
69 void error(
double e)
const { last_error_ = e; }
78 mutable double last_error_;
91 template <
typename MatrixType,
typename VectorType>
96 VectorType result = rhs;
99 VectorType residual = rhs;
101 VectorType r0star = rhs;
102 VectorType tmp0 = rhs;
103 VectorType tmp1 = rhs;
107 CPU_ScalarType ip_rr0star = norm_rhs_host * norm_rhs_host;
109 CPU_ScalarType alpha;
110 CPU_ScalarType omega;
112 CPU_ScalarType new_ip_rr0star = 0;
113 CPU_ScalarType residual_norm = norm_rhs_host;
115 if (norm_rhs_host == 0)
118 bool restart_flag =
true;
129 ip_rr0star *= ip_rr0star;
130 restart_flag =
false;
138 s = residual - alpha*tmp0;
144 result += alpha * p + omega * s;
145 residual = s - omega * tmp1;
149 if (std::fabs(residual_norm / norm_rhs_host) < tag.
tolerance())
152 beta = new_ip_rr0star / ip_rr0star * alpha/omega;
153 ip_rr0star = new_ip_rr0star;
162 p = residual + beta * p;
166 tag.
error(residual_norm / norm_rhs_host);
171 template <
typename MatrixType,
typename VectorType>
174 return solve(matrix, rhs, tag);
187 template <
typename MatrixType,
typename VectorType,
typename PreconditionerType>
188 VectorType
solve(
const MatrixType &
matrix, VectorType
const & rhs,
bicgstab_tag const & tag, PreconditionerType
const & precond)
192 VectorType result = rhs;
195 VectorType residual = rhs;
196 VectorType r0star = residual;
197 VectorType tmp0 = rhs;
198 VectorType tmp1 = rhs;
201 VectorType p = residual;
206 CPU_ScalarType alpha;
207 CPU_ScalarType omega;
208 CPU_ScalarType new_ip_rr0star = 0;
209 CPU_ScalarType residual_norm = norm_rhs_host;
211 if (norm_rhs_host == 0)
214 bool restart_flag =
true;
222 precond.apply(residual);
226 ip_rr0star *= ip_rr0star;
227 restart_flag =
false;
236 s = residual - alpha*tmp0;
243 result += alpha * p + omega * s;
244 residual = s - omega * tmp1;
247 if (residual_norm / norm_rhs_host < tag.
tolerance())
252 beta = new_ip_rr0star / ip_rr0star * alpha/omega;
253 ip_rr0star = new_ip_rr0star;
262 p = residual + beta * p;
268 tag.
error(residual_norm / norm_rhs_host);
std::size_t vcl_size_t
Definition: forwards.h:58
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:86
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
Generic size and resize functionality for different vector and matrix types.
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
A dense matrix class.
Definition: forwards.h:293
void error(double e) const
Sets the estimated relative error at the end of the solver run.
Definition: bicgstab.hpp:69
void clear(VectorType &vec)
Generic routine for setting all entries of a vector to zero. This is the version for non-ViennaCL obj...
Definition: clear.hpp:57
This file provides the forward declarations for the main types used within ViennaCL.
vcl_size_t max_iterations_before_restart() const
Returns the maximum number of iterations before a restart.
Definition: bicgstab.hpp:60
viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< VectorT1 >::type >::value, typename VectorT1::value_type >::type inner_prod(VectorT1 const &v1, VectorT2 const &v2)
Definition: inner_prod.hpp:89
Generic interface for the computation of inner products. See viennacl/linalg/vector_operations.hpp for implementations.
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
bicgstab_tag(double tol=1e-8, vcl_size_t max_iters=400, vcl_size_t max_iters_before_restart=200)
The constructor.
Definition: bicgstab.hpp:52
A tag class representing the use of no preconditioner.
Definition: forwards.h:727
void iters(vcl_size_t i) const
Definition: bicgstab.hpp:64
T::value_type type
Definition: result_of.hpp:230
Generic clear functionality for different vector and matrix types.
double error() const
Returns the estimated relative error at the end of the solver run.
Definition: bicgstab.hpp:67
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
A tag for the stabilized Bi-conjugate gradient solver. Used for supplying solver parameters and for d...
Definition: bicgstab.hpp:43
vcl_size_t iters() const
Return the number of solver iterations:
Definition: bicgstab.hpp:63
VectorType solve(const MatrixType &matrix, VectorType const &rhs, bicgstab_tag const &tag)
Implementation of the stabilized Bi-conjugate gradient solver.
Definition: bicgstab.hpp:92
A collection of compile time type deductions.
vcl_size_t max_iterations() const
Returns the maximum number of iterations.
Definition: bicgstab.hpp:58
double tolerance() const
Returns the relative tolerance.
Definition: bicgstab.hpp:56