ViennaCL - The Vienna Computing Library  1.5.1
cg.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_CG_HPP_
2 #define VIENNACL_LINALG_CG_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 <map>
27 #include <cmath>
28 #include "viennacl/forwards.h"
29 #include "viennacl/tools/tools.hpp"
30 #include "viennacl/linalg/ilu.hpp"
31 #include "viennacl/linalg/prod.hpp"
35 #include "viennacl/traits/size.hpp"
37 
38 namespace viennacl
39 {
40  namespace linalg
41  {
42 
45  class cg_tag
46  {
47  public:
53  cg_tag(double tol = 1e-8, unsigned int max_iterations = 300) : tol_(tol), iterations_(max_iterations) {}
54 
56  double tolerance() const { return tol_; }
58  unsigned int max_iterations() const { return iterations_; }
59 
61  unsigned int iters() const { return iters_taken_; }
62  void iters(unsigned int i) const { iters_taken_ = i; }
63 
65  double error() const { return last_error_; }
67  void error(double e) const { last_error_ = e; }
68 
69 
70  private:
71  double tol_;
72  unsigned int iterations_;
73 
74  //return values from solver
75  mutable unsigned int iters_taken_;
76  mutable double last_error_;
77  };
78 
79 
89  template <typename MatrixType, typename VectorType>
90  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag)
91  {
92  //typedef typename VectorType::value_type ScalarType;
93  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
94  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
95  //std::cout << "Starting CG" << std::endl;
96  VectorType result = rhs;
98 
99  VectorType residual = rhs;
100  VectorType p = rhs;
101  VectorType tmp = rhs;
102 
103  CPU_ScalarType ip_rr = viennacl::linalg::inner_prod(rhs,rhs);
104  CPU_ScalarType alpha;
105  CPU_ScalarType new_ip_rr = 0;
106  CPU_ScalarType beta;
107  CPU_ScalarType norm_rhs = std::sqrt(ip_rr);
108 
109  //std::cout << "Starting CG solver iterations... " << std::endl;
110  if (norm_rhs == 0) //solution is zero if RHS norm is zero
111  return result;
112 
113  for (unsigned int i = 0; i < tag.max_iterations(); ++i)
114  {
115  tag.iters(i+1);
116  tmp = viennacl::linalg::prod(matrix, p);
117 
118  alpha = ip_rr / viennacl::linalg::inner_prod(tmp, p);
119  result += alpha * p;
120  residual -= alpha * tmp;
121 
122  new_ip_rr = viennacl::linalg::norm_2(residual);
123  if (new_ip_rr / norm_rhs < tag.tolerance())
124  break;
125  new_ip_rr *= new_ip_rr;
126 
127  beta = new_ip_rr / ip_rr;
128  ip_rr = new_ip_rr;
129 
130  p = residual + beta * p;
131  }
132 
133  //store last error estimate:
134  tag.error(std::sqrt(new_ip_rr) / norm_rhs);
135 
136  return result;
137  }
138 
139  template <typename MatrixType, typename VectorType>
140  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag, viennacl::linalg::no_precond)
141  {
142  return solve(matrix, rhs, tag);
143  }
144 
155  template <typename MatrixType, typename VectorType, typename PreconditionerType>
156  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag, PreconditionerType const & precond)
157  {
158  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
159  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
160 
161  VectorType result = rhs;
162  viennacl::traits::clear(result);
163 
164  VectorType residual = rhs;
165  VectorType tmp = rhs;
166  VectorType z = rhs;
167 
168  precond.apply(z);
169  VectorType p = z;
170 
171  CPU_ScalarType ip_rr = viennacl::linalg::inner_prod(residual, z);
172  CPU_ScalarType alpha;
173  CPU_ScalarType new_ip_rr = 0;
174  CPU_ScalarType beta;
175  CPU_ScalarType norm_rhs_squared = ip_rr;
176  CPU_ScalarType new_ipp_rr_over_norm_rhs;
177 
178  if (norm_rhs_squared == 0) //solution is zero if RHS norm is zero
179  return result;
180 
181  for (unsigned int i = 0; i < tag.max_iterations(); ++i)
182  {
183  tag.iters(i+1);
184  tmp = viennacl::linalg::prod(matrix, p);
185 
186  alpha = ip_rr / viennacl::linalg::inner_prod(tmp, p);
187 
188  result += alpha * p;
189  residual -= alpha * tmp;
190  z = residual;
191  precond.apply(z);
192 
193  new_ip_rr = viennacl::linalg::inner_prod(residual, z);
194  new_ipp_rr_over_norm_rhs = new_ip_rr / norm_rhs_squared;
195  if (std::fabs(new_ipp_rr_over_norm_rhs) < tag.tolerance() * tag.tolerance()) //squared norms involved here
196  break;
197 
198  beta = new_ip_rr / ip_rr;
199  ip_rr = new_ip_rr;
200 
201  p = z + beta*p;
202  }
203 
204  //store last error estimate:
205  tag.error(std::sqrt(std::fabs(new_ip_rr / norm_rhs_squared)));
206 
207  return result;
208  }
209 
210  }
211 }
212 
213 #endif
T norm_2(std::vector< T, A > const &v1)
Definition: norm_2.hpp:86
void error(double e) const
Sets the estimated relative error at the end of the solver run.
Definition: cg.hpp:67
double tolerance() const
Returns the relative tolerance.
Definition: cg.hpp:56
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...
cg_tag(double tol=1e-8, unsigned int max_iterations=300)
The constructor.
Definition: cg.hpp:53
Various little tools used here and there in ViennaCL.
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
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
A tag class representing the use of no preconditioner.
Definition: forwards.h:727
Implementations of incomplete factorization preconditioners. Convenience header file.
void iters(unsigned int i) const
Definition: cg.hpp:62
unsigned int iters() const
Return the number of solver iterations:
Definition: cg.hpp:61
T::value_type type
Definition: result_of.hpp:230
Generic clear functionality for different vector and matrix types.
unsigned int max_iterations() const
Returns the maximum number of iterations.
Definition: cg.hpp:58
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
A tag for the conjugate gradient Used for supplying solver parameters and for dispatching the solve()...
Definition: cg.hpp:45
double error() const
Returns the estimated relative error at the end of the solver run.
Definition: cg.hpp:65
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.