ViennaCL - The Vienna Computing Library  1.2.0
cg.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_CG_HPP_
2 #define VIENNACL_CG_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
24 #include <vector>
25 #include <map>
26 #include <cmath>
27 #include "viennacl/forwards.h"
28 #include "viennacl/tools/tools.hpp"
29 #include "viennacl/linalg/ilu.hpp"
30 #include "viennacl/linalg/prod.hpp"
33 #include "viennacl/traits/size.hpp"
35 
36 namespace viennacl
37 {
38  namespace linalg
39  {
40 
43  class cg_tag
44  {
45  public:
51  cg_tag(double tol = 1e-8, unsigned int max_iterations = 300) : _tol(tol), _iterations(max_iterations) {};
52 
54  double tolerance() const { return _tol; }
56  unsigned int max_iterations() const { return _iterations; }
57 
59  unsigned int iters() const { return iters_taken_; }
60  void iters(unsigned int i) const { iters_taken_ = i; }
61 
63  double error() const { return last_error_; }
65  void error(double e) const { last_error_ = e; }
66 
67 
68  private:
69  double _tol;
70  unsigned int _iterations;
71 
72  //return values from solver
73  mutable unsigned int iters_taken_;
74  mutable double last_error_;
75  };
76 
77 
87  template <typename MatrixType, typename VectorType>
88  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag)
89  {
90  //typedef typename VectorType::value_type ScalarType;
91  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
92  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
93  //std::cout << "Starting CG" << std::endl;
94  std::size_t problem_size = viennacl::traits::size(rhs);
95  VectorType result(problem_size);
97 
98  VectorType residual = rhs;
99  VectorType p = rhs;
100  VectorType tmp(problem_size);
101 
102  ScalarType tmp_in_p;
103  ScalarType residual_norm_squared;
104  CPU_ScalarType ip_rr = viennacl::linalg::inner_prod(rhs,rhs);
105  CPU_ScalarType alpha;
106  CPU_ScalarType new_ip_rr = 0;
107  CPU_ScalarType beta;
108  CPU_ScalarType norm_rhs_squared = ip_rr;
109 
110  //std::cout << "Starting CG solver iterations... " << std::endl;
111 
112  for (unsigned int i = 0; i < tag.max_iterations(); ++i)
113  {
114  tag.iters(i+1);
115  tmp = viennacl::linalg::prod(matrix, p);
116 
117  //alpha = ip_rr / viennacl::linalg::inner_prod(tmp, p);
118  tmp_in_p = viennacl::linalg::inner_prod(tmp, p);
119  alpha = ip_rr / static_cast<CPU_ScalarType>(tmp_in_p);
120 
121  result += alpha * p;
122  residual -= alpha * tmp;
123 
124  residual_norm_squared = viennacl::linalg::inner_prod(residual,residual);
125  new_ip_rr = static_cast<CPU_ScalarType>(residual_norm_squared);
126  if (new_ip_rr / norm_rhs_squared < tag.tolerance() * tag.tolerance())//squared norms involved here
127  break;
128 
129  beta = new_ip_rr / ip_rr;
130  ip_rr = new_ip_rr;
131 
132  //p = residual + beta*p;
133  p *= beta;
134  p += residual;
135  }
136 
137  //store last error estimate:
138  tag.error(sqrt(new_ip_rr / norm_rhs_squared));
139 
140  return result;
141  }
142 
143  template <typename MatrixType, typename VectorType>
144  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag, viennacl::linalg::no_precond)
145  {
146  return solve(matrix, rhs, tag);
147  }
148 
159  template <typename MatrixType, typename VectorType, typename PreconditionerType>
160  VectorType solve(const MatrixType & matrix, VectorType const & rhs, cg_tag const & tag, PreconditionerType const & precond)
161  {
162  typedef typename viennacl::result_of::value_type<VectorType>::type ScalarType;
163  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
164  unsigned int problem_size = viennacl::traits::size(rhs);
165 
166  VectorType result(problem_size);
167  result.clear();
168 
169  VectorType residual = rhs;
170  VectorType tmp(problem_size);
171  VectorType z = rhs;
172 
173  precond.apply(z);
174  VectorType p = z;
175 
176  ScalarType tmp_in_p;
177  ScalarType residual_in_z;
178  CPU_ScalarType ip_rr = viennacl::linalg::inner_prod(residual, z);
179  CPU_ScalarType alpha;
180  CPU_ScalarType new_ip_rr = 0;
181  CPU_ScalarType beta;
182  CPU_ScalarType norm_rhs_squared = ip_rr;
183  CPU_ScalarType new_ipp_rr_over_norm_rhs;
184 
185  for (unsigned int i = 0; i < tag.max_iterations(); ++i)
186  {
187  tag.iters(i+1);
188  tmp = viennacl::linalg::prod(matrix, p);
189 
190  tmp_in_p = viennacl::linalg::inner_prod(tmp, p);
191  alpha = ip_rr / static_cast<CPU_ScalarType>(tmp_in_p);
192 
193  result += alpha * p;
194  residual -= alpha * tmp;
195  z = residual;
196  precond.apply(z);
197 
198  residual_in_z = viennacl::linalg::inner_prod(residual, z);
199  new_ip_rr = static_cast<CPU_ScalarType>(residual_in_z);
200  new_ipp_rr_over_norm_rhs = new_ip_rr / norm_rhs_squared;
201  if (std::fabs(new_ipp_rr_over_norm_rhs) < tag.tolerance() * tag.tolerance()) //squared norms involved here
202  break;
203 
204  beta = new_ip_rr / ip_rr;
205  ip_rr = new_ip_rr;
206 
207  //p = z + beta*p;
208  p *= beta;
209  p += z;
210  }
211 
212  //store last error estimate:
213  tag.error(sqrt(std::fabs(new_ip_rr / norm_rhs_squared)));
214 
215  return result;
216  }
217 
218  }
219 }
220 
221 #endif