ViennaCL - The Vienna Computing Library  1.5.1
bicgstab.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_BICGSTAB_HPP_
2 #define VIENNACL_LINALG_BICGSTAB_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
25 #include <vector>
26 #include <cmath>
27 #include "viennacl/forwards.h"
28 #include "viennacl/tools/tools.hpp"
29 #include "viennacl/linalg/prod.hpp"
33 #include "viennacl/traits/size.hpp"
35 
36 namespace viennacl
37 {
38  namespace linalg
39  {
40 
44  {
45  public:
52  bicgstab_tag(double tol = 1e-8, vcl_size_t max_iters = 400, vcl_size_t max_iters_before_restart = 200)
53  : tol_(tol), iterations_(max_iters), iterations_before_restart_(max_iters_before_restart) {}
54 
56  double tolerance() const { return tol_; }
58  vcl_size_t max_iterations() const { return iterations_; }
60  vcl_size_t max_iterations_before_restart() const { return iterations_before_restart_; }
61 
63  vcl_size_t iters() const { return iters_taken_; }
64  void iters(vcl_size_t i) const { iters_taken_ = i; }
65 
67  double error() const { return last_error_; }
69  void error(double e) const { last_error_ = e; }
70 
71  private:
72  double tol_;
73  vcl_size_t iterations_;
74  vcl_size_t iterations_before_restart_;
75 
76  //return values from solver
77  mutable vcl_size_t iters_taken_;
78  mutable double last_error_;
79  };
80 
81 
91  template <typename MatrixType, typename VectorType>
92  VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag)
93  {
94  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
95  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
96  VectorType result = rhs;
98 
99  VectorType residual = rhs;
100  VectorType p = rhs;
101  VectorType r0star = rhs;
102  VectorType tmp0 = rhs;
103  VectorType tmp1 = rhs;
104  VectorType s = rhs;
105 
106  CPU_ScalarType norm_rhs_host = viennacl::linalg::norm_2(residual);
107  CPU_ScalarType ip_rr0star = norm_rhs_host * norm_rhs_host;
108  CPU_ScalarType beta;
109  CPU_ScalarType alpha;
110  CPU_ScalarType omega;
111  //ScalarType inner_prod_temp; //temporary variable for inner product computation
112  CPU_ScalarType new_ip_rr0star = 0;
113  CPU_ScalarType residual_norm = norm_rhs_host;
114 
115  if (norm_rhs_host == 0) //solution is zero if RHS norm is zero
116  return result;
117 
118  bool restart_flag = true;
119  vcl_size_t last_restart = 0;
120  for (vcl_size_t i = 0; i < tag.max_iterations(); ++i)
121  {
122  if (restart_flag)
123  {
124  residual = rhs;
125  residual -= viennacl::linalg::prod(matrix, result);
126  p = residual;
127  r0star = residual;
128  ip_rr0star = viennacl::linalg::norm_2(residual);
129  ip_rr0star *= ip_rr0star;
130  restart_flag = false;
131  last_restart = i;
132  }
133 
134  tag.iters(i+1);
135  tmp0 = viennacl::linalg::prod(matrix, p);
136  alpha = ip_rr0star / viennacl::linalg::inner_prod(tmp0, r0star);
137 
138  s = residual - alpha*tmp0;
139 
140  tmp1 = viennacl::linalg::prod(matrix, s);
141  CPU_ScalarType norm_tmp1 = viennacl::linalg::norm_2(tmp1);
142  omega = viennacl::linalg::inner_prod(tmp1, s) / (norm_tmp1 * norm_tmp1);
143 
144  result += alpha * p + omega * s;
145  residual = s - omega * tmp1;
146 
147  new_ip_rr0star = viennacl::linalg::inner_prod(residual, r0star);
148  residual_norm = viennacl::linalg::norm_2(residual);
149  if (std::fabs(residual_norm / norm_rhs_host) < tag.tolerance())
150  break;
151 
152  beta = new_ip_rr0star / ip_rr0star * alpha/omega;
153  ip_rr0star = new_ip_rr0star;
154 
155  if (ip_rr0star == 0 || omega == 0 || i - last_restart > tag.max_iterations_before_restart()) //search direction degenerate. A restart might help
156  restart_flag = true;
157 
158  // Execution of
159  // p = residual + beta * (p - omega*tmp0);
160  // without introducing temporary vectors:
161  p -= omega * tmp0;
162  p = residual + beta * p;
163  }
164 
165  //store last error estimate:
166  tag.error(residual_norm / norm_rhs_host);
167 
168  return result;
169  }
170 
171  template <typename MatrixType, typename VectorType>
172  VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag, viennacl::linalg::no_precond)
173  {
174  return solve(matrix, rhs, tag);
175  }
176 
187  template <typename MatrixType, typename VectorType, typename PreconditionerType>
188  VectorType solve(const MatrixType & matrix, VectorType const & rhs, bicgstab_tag const & tag, PreconditionerType const & precond)
189  {
190  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
191  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
192  VectorType result = rhs;
193  viennacl::traits::clear(result);
194 
195  VectorType residual = rhs;
196  VectorType r0star = residual; //can be chosen arbitrarily in fact
197  VectorType tmp0 = rhs;
198  VectorType tmp1 = rhs;
199  VectorType s = rhs;
200 
201  VectorType p = residual;
202 
203  CPU_ScalarType ip_rr0star = viennacl::linalg::norm_2(residual);
204  CPU_ScalarType norm_rhs_host = viennacl::linalg::norm_2(residual);
205  CPU_ScalarType beta;
206  CPU_ScalarType alpha;
207  CPU_ScalarType omega;
208  CPU_ScalarType new_ip_rr0star = 0;
209  CPU_ScalarType residual_norm = norm_rhs_host;
210 
211  if (norm_rhs_host == 0) //solution is zero if RHS norm is zero
212  return result;
213 
214  bool restart_flag = true;
215  vcl_size_t last_restart = 0;
216  for (unsigned int i = 0; i < tag.max_iterations(); ++i)
217  {
218  if (restart_flag)
219  {
220  residual = rhs;
221  residual -= viennacl::linalg::prod(matrix, result);
222  precond.apply(residual);
223  p = residual;
224  r0star = residual;
225  ip_rr0star = viennacl::linalg::norm_2(residual);
226  ip_rr0star *= ip_rr0star;
227  restart_flag = false;
228  last_restart = i;
229  }
230 
231  tag.iters(i+1);
232  tmp0 = viennacl::linalg::prod(matrix, p);
233  precond.apply(tmp0);
234  alpha = ip_rr0star / viennacl::linalg::inner_prod(tmp0, r0star);
235 
236  s = residual - alpha*tmp0;
237 
238  tmp1 = viennacl::linalg::prod(matrix, s);
239  precond.apply(tmp1);
240  CPU_ScalarType norm_tmp1 = viennacl::linalg::norm_2(tmp1);
241  omega = viennacl::linalg::inner_prod(tmp1, s) / (norm_tmp1 * norm_tmp1);
242 
243  result += alpha * p + omega * s;
244  residual = s - omega * tmp1;
245 
246  residual_norm = viennacl::linalg::norm_2(residual);
247  if (residual_norm / norm_rhs_host < tag.tolerance())
248  break;
249 
250  new_ip_rr0star = viennacl::linalg::inner_prod(residual, r0star);
251 
252  beta = new_ip_rr0star / ip_rr0star * alpha/omega;
253  ip_rr0star = new_ip_rr0star;
254 
255  if (ip_rr0star == 0 || omega == 0 || i - last_restart > tag.max_iterations_before_restart()) //search direction degenerate. A restart might help
256  restart_flag = true;
257 
258  // Execution of
259  // p = residual + beta * (p - omega*tmp0);
260  // without introducing temporary vectors:
261  p -= omega * tmp0;
262  p = residual + beta * p;
263 
264  //std::cout << "Rel. Residual in current step: " << std::sqrt(std::fabs(viennacl::linalg::inner_prod(residual, residual) / norm_rhs_host)) << std::endl;
265  }
266 
267  //store last error estimate:
268  tag.error(residual_norm / norm_rhs_host);
269 
270  return result;
271  }
272 
273  }
274 }
275 
276 #endif
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...
Various little tools used here and there in ViennaCL.
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.
A dense matrix class.
Definition: forwards.h:293
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.
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