ViennaCL - The Vienna Computing Library  1.2.0
row_scaling.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_ROW_SCALING_HPP_
2 #define VIENNACL_LINALG_ROW_SCALING_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 
41  {
42  public:
47  row_scaling_tag(unsigned int p = 2) : norm_(p) {}
48 
50  unsigned int norm() const { return norm_; }
51 
52  private:
53  unsigned int norm_;
54  };
55 
56 
59  template <typename MatrixType>
61  {
62  typedef typename MatrixType::value_type ScalarType;
63 
64  public:
70  row_scaling(MatrixType const & mat, row_scaling_tag const & tag) : system_matrix(mat), tag_(tag)
71  {
72  assert(mat.size1() == mat.size2());
73  diag_M_inv.resize(mat.size1()); //resize without preserving values
74 
75  for (typename MatrixType::const_iterator1 row_it = system_matrix.begin1();
76  row_it != system_matrix.end1();
77  ++row_it)
78  {
79  diag_M_inv[row_it.index1()];
80  for (typename MatrixType::const_iterator2 col_it = row_it.begin();
81  col_it != row_it.end();
82  ++col_it)
83  {
84  if (tag_.norm() == 1)
85  diag_M_inv[col_it.index1()] += std::fabs(*col_it);
86  else
87  diag_M_inv[col_it.index1()] += (*col_it) * (*col_it);
88  }
89  if (diag_M_inv[row_it.index1()] == 0)
90  throw "ViennaCL: Zero row encountered while setting up row scaling preconditioner!";
91 
92  if (tag_.norm() == 1)
93  diag_M_inv[row_it.index1()] = static_cast<ScalarType>(1.0) / diag_M_inv[row_it.index1()];
94  else
95  diag_M_inv[row_it.index1()] = static_cast<ScalarType>(1.0) / std::sqrt(diag_M_inv[row_it.index1()]);
96  }
97  }
98 
99 
101  template <typename VectorType>
102  void apply(VectorType & vec) const
103  {
104  assert(vec.size() == diag_M_inv.size());
105  for (size_t i=0; i<vec.size(); ++i)
106  {
107  vec[i] *= diag_M_inv[i];
108  }
109  }
110 
111  private:
112  MatrixType const & system_matrix;
113  row_scaling_tag const & tag_;
114  std::vector<ScalarType> diag_M_inv;
115  };
116 
117 
122  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
123  class row_scaling< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
124  {
126 
127  public:
133  row_scaling(MatrixType const & mat, row_scaling_tag const & tag) : system_matrix(mat), tag_(tag), diag_M_inv(mat.size1())
134  {
135  assert(system_matrix.size1() == system_matrix.size2());
136 
137  init_gpu();
138  }
139 
140  /*
141  void init_cpu()
142  {
143  std::vector< std::map<unsigned int, ScalarType> > cpu_check;
144  std::vector<ScalarType> diag_M_inv_cpu(system_matrix.size1());
145 
146  copy(system_matrix, cpu_check);
147  viennacl::tools::const_sparse_matrix_adapter<ScalarType> cpu_check_adapter(cpu_check);
148 
149  for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator1 row_it = cpu_check_adapter.begin1();
150  row_it != cpu_check_adapter.end1();
151  ++row_it)
152  {
153  diag_M_inv_cpu[row_it.index1()] = 0;
154  for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator2 col_it = row_it.begin();
155  col_it != row_it.end();
156  ++col_it)
157  {
158  if (tag_.norm() == 1)
159  diag_M_inv_cpu[col_it.index1()] += std::fabs(*col_it);
160  else
161  diag_M_inv_cpu[col_it.index1()] += (*col_it) * (*col_it);
162  }
163  if (diag_M_inv_cpu[row_it.index1()] == 0)
164  throw "ViennaCL: Zero row encountered while setting up row scaling preconditioner!";
165 
166  if (tag_.norm() == 1)
167  diag_M_inv_cpu[row_it.index1()] = static_cast<ScalarType>(1.0) / diag_M_inv_cpu[row_it.index1()];
168  else
169  diag_M_inv_cpu[row_it.index1()] = static_cast<ScalarType>(1.0) / std::sqrt(diag_M_inv_cpu[row_it.index1()]);
170  }
171 
172  diag_M_inv.resize(system_matrix.size1(), false);
173  viennacl::fast_copy(diag_M_inv_cpu, diag_M_inv);
174  } */
175 
176  void init_gpu()
177  {
178  if (tag_.norm() == 1)
179  {
181  viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
182  "row_scaling_1");
183 
184  viennacl::ocl::enqueue( k(system_matrix.handle1(), system_matrix.handle2(), system_matrix.handle(),
185  diag_M_inv, static_cast<cl_uint>(diag_M_inv.size())) );
186  }
187  else
188  {
190  viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
191  "row_scaling_2");
192 
193  viennacl::ocl::enqueue( k(system_matrix.handle1(), system_matrix.handle2(), system_matrix.handle(),
194  diag_M_inv, static_cast<cl_uint>(diag_M_inv.size())) );
195  }
196  }
197 
198  template <unsigned int ALIGNMENT>
200  {
201  assert(viennacl::traits::size1(system_matrix) == viennacl::traits::size(vec));
202 
203  //run kernel:
204  viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::vector<ScalarType, ALIGNMENT>::program_name(),
205  "diag_precond");
206 
208  k(viennacl::traits::handle(diag_M_inv), cl_uint(viennacl::traits::start(diag_M_inv)), cl_uint(viennacl::traits::size(diag_M_inv)),
209  viennacl::traits::handle(vec), cl_uint(viennacl::traits::start(vec)), cl_uint(viennacl::traits::size(vec)) )
210  );
211 
212  }
213 
214  private:
215  MatrixType const & system_matrix;
216  row_scaling_tag const & tag_;
217  viennacl::vector<ScalarType> diag_M_inv;
218  };
219 
220  }
221 }
222 
223 
224 
225 
226 #endif
227 
228 
229