ViennaCL - The Vienna Computing Library  1.2.0
jacobi_precond.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_JACOBI_PRECOND_HPP_
2 #define VIENNACL_JACOBI_PRECOND_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 <cmath>
26 #include "viennacl/forwards.h"
27 #include "viennacl/vector.hpp"
29 #include "viennacl/tools/tools.hpp"
30 
31 #include <map>
32 
33 namespace viennacl
34 {
35  namespace linalg
36  {
37 
40  class jacobi_tag {};
41 
42 
45  template <typename MatrixType>
47  {
48  typedef typename MatrixType::value_type ScalarType;
49 
50  public:
51  jacobi_precond(MatrixType const & mat, jacobi_tag const & tag) : system_matrix(mat)
52  {
53  assert(mat.size1() == mat.size2());
54  diag_A_inv.resize(mat.size1()); //resize without preserving values
55 
56  for (typename MatrixType::const_iterator1 row_it = system_matrix.begin1();
57  row_it != system_matrix.end1();
58  ++row_it)
59  {
60  bool diag_found = false;
61  for (typename MatrixType::const_iterator2 col_it = row_it.begin();
62  col_it != row_it.end();
63  ++col_it)
64  {
65  if (col_it.index1() == col_it.index2())
66  {
67  diag_A_inv[col_it.index1()] = static_cast<ScalarType>(1.0) / *col_it;
68  diag_found = true;
69  }
70  }
71  if (!diag_found)
72  throw "ViennaCL: Zero in diagonal encountered while setting up Jacobi preconditioner!";
73  }
74  }
75 
76 
78  template <typename VectorType>
79  void apply(VectorType & vec) const
80  {
81  assert(vec.size() == diag_A_inv.size());
82  for (size_t i=0; i<vec.size(); ++i)
83  {
84  vec[i] *= diag_A_inv[i];
85  }
86  }
87 
88  private:
89  MatrixType const & system_matrix;
90  std::vector<ScalarType> diag_A_inv;
91  };
92 
93 
98  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
99  class jacobi_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
100  {
102 
103  public:
104  jacobi_precond(MatrixType const & mat, jacobi_tag const & tag) : system_matrix(mat), diag_A_inv(mat.size1())
105  {
106  assert(system_matrix.size1() == system_matrix.size2());
107 
108  init_gpu();
109  }
110 
111  /*void init_cpu()
112  {
113 
114  std::vector< std::map<unsigned int, ScalarType> > cpu_check;
115  std::vector<ScalarType> diag_A_inv_cpu(system_matrix.size1());
116 
117  copy(system_matrix, cpu_check);
118  viennacl::tools::const_sparse_matrix_adapter<ScalarType> cpu_check_adapter(cpu_check);
119 
120  for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator1 row_it = cpu_check_adapter.begin1();
121  row_it != cpu_check_adapter.end1();
122  ++row_it)
123  {
124  bool diag_found = false;
125  for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator2 col_it = row_it.begin();
126  col_it != row_it.end();
127  ++col_it)
128  {
129  if (col_it.index1() == col_it.index2())
130  {
131  diag_found = true;
132  diag_A_inv_cpu[col_it.index1()] = static_cast<ScalarType>(1.0) / *col_it;
133  }
134  }
135  if (!diag_found)
136  throw "ViennaCL: Zero in diagonal encountered while setting up Jacobi preconditioner!";
137  }
138 
139  diag_A_inv.resize(system_matrix.size1(), false);
140  viennacl::fast_copy(diag_A_inv_cpu, diag_A_inv);
141  }*/
142 
143  void init_gpu()
144  {
146  viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
147  "jacobi_precond");
148 
149  viennacl::ocl::enqueue( k(system_matrix.handle1(), system_matrix.handle2(), system_matrix.handle(),
150  diag_A_inv, static_cast<cl_uint>(diag_A_inv.size())) );
151  }
152 
153 
154  template <unsigned int ALIGNMENT>
156  {
157  assert(viennacl::traits::size1(system_matrix) == viennacl::traits::size(vec));
158 
159  //run kernel:
160  viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::vector<ScalarType, ALIGNMENT>::program_name(),
161  "diag_precond");
162 
164  k(viennacl::traits::handle(diag_A_inv), cl_uint(viennacl::traits::start(diag_A_inv)), cl_uint(viennacl::traits::size(diag_A_inv)),
165  viennacl::traits::handle(vec), cl_uint(viennacl::traits::start(vec)), cl_uint(viennacl::traits::size(vec)) )
166  );
167  }
168 
169  private:
170  MatrixType const & system_matrix;
171  viennacl::vector<ScalarType> diag_A_inv;
172  };
173 
174  }
175 }
176 
177 
178 
179 
180 #endif
181 
182 
183