ViennaCL - The Vienna Computing Library  1.5.1
qr-method-common.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_QR_METHOD_COMMON_HPP
2 #define VIENNACL_LINALG_QR_METHOD_COMMON_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 
21 #include <cmath>
22 
23 #include "viennacl/ocl/device.hpp"
24 #include "viennacl/ocl/handle.hpp"
25 #include "viennacl/ocl/kernel.hpp"
28 #include "viennacl/vector.hpp"
29 #include "viennacl/matrix.hpp"
30 
31 #include <boost/numeric/ublas/vector.hpp>
32 #include <boost/numeric/ublas/io.hpp>
33 
38 namespace viennacl
39 {
40  namespace linalg
41  {
42  const std::string SVD_BIDIAG_PACK_KERNEL = "bidiag_pack";
43  const std::string SVD_HOUSEHOLDER_UPDATE_A_LEFT_KERNEL = "house_update_A_left";
44  const std::string SVD_HOUSEHOLDER_UPDATE_A_RIGHT_KERNEL = "house_update_A_right";
45  const std::string SVD_HOUSEHOLDER_UPDATE_QL_KERNEL = "house_update_QL";
46  const std::string SVD_HOUSEHOLDER_UPDATE_QR_KERNEL = "house_update_QR";
47  const std::string SVD_COPY_COL_KERNEL = "copy_col";
48  const std::string SVD_COPY_ROW_KERNEL = "copy_row";
49  const std::string SVD_MATRIX_TRANSPOSE_KERNEL = "transpose_inplace";
50  const std::string SVD_INVERSE_SIGNS_KERNEL = "inverse_signs";
51  const std::string SVD_GIVENS_PREV_KERNEL = "givens_prev";
52  const std::string SVD_GIVENS_NEXT_KERNEL = "givens_next";
53  const std::string SVD_FINAL_ITER_UPDATE_KERNEL = "final_iter_update";
54  const std::string SVD_UPDATE_QR_COLUMN_KERNEL = "update_qr_column";
55 
56  namespace detail
57  {
58  //static const float EPS = 0.00001f;
59  //static const vcl_size_t ITER_MAX = 50;
60 
61  static const double EPS = 1e-10;
62  static const vcl_size_t ITER_MAX = 50;
63 
64  template <typename SCALARTYPE>
65  SCALARTYPE pythag(SCALARTYPE a, SCALARTYPE b)
66  {
67  return std::sqrt(a*a + b*b);
68  }
69 
70  template <typename SCALARTYPE>
71  SCALARTYPE sign(SCALARTYPE val)
72  {
73  return (val >= 0) ? SCALARTYPE(1) : SCALARTYPE(-1);
74  }
75 
76  // DEPRECATED: Replace with viennacl::linalg::norm_2
77  template <typename VectorType>
78  typename VectorType::value_type norm_lcl(VectorType const & x, vcl_size_t size)
79  {
80  typename VectorType::value_type x_norm = 0.0;
81  for(vcl_size_t i = 0; i < size; i++)
82  x_norm += std::pow(x[i], 2);
83  return std::sqrt(x_norm);
84  }
85 
86  template <typename VectorType>
87  void normalize(VectorType & x, vcl_size_t size)
88  {
89  typename VectorType::value_type x_norm = norm_lcl(x, size);
90  for(vcl_size_t i = 0; i < size; i++)
91  x[i] /= x_norm;
92  }
93 
94 
95 
96  template <typename VectorType>
97  void householder_vector(VectorType & v, vcl_size_t start)
98  {
99  typedef typename VectorType::value_type ScalarType;
100  ScalarType x_norm = norm_lcl(v, v.size());
101  ScalarType alpha = -sign(v[start]) * x_norm;
102  v[start] += alpha;
103  normalize(v, v.size());
104  }
105 
106  template <typename MatrixType>
107  void transpose(MatrixType & A)
108  {
109  typedef typename MatrixType::value_type ScalarType;
110  typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
111 
113 
114  viennacl::ocl::enqueue(kernel(A,
115  static_cast<cl_uint>(A.internal_size1()),
116  static_cast<cl_uint>(A.internal_size2())
117  )
118  );
119  }
120 
121 
122 
123  template <typename T>
124  void cdiv(T xr, T xi, T yr, T yi, T& cdivr, T& cdivi)
125  {
126  // Complex scalar division.
127  T r;
128  T d;
129  if (std::fabs(yr) > std::fabs(yi))
130  {
131  r = yi / yr;
132  d = yr + r * yi;
133  cdivr = (xr + r * xi) / d;
134  cdivi = (xi - r * xr) / d;
135  }
136  else
137  {
138  r = yr / yi;
139  d = yi + r * yr;
140  cdivr = (r * xr + xi) / d;
141  cdivi = (r * xi - xr) / d;
142  }
143  }
144 
145 
146  template <typename SCALARTYPE, unsigned int ALIGNMENT>
149  vcl_size_t row_start,
150  vcl_size_t col_start,
151  bool copy_col
152  )
153  {
154 
155  std::string kernel_name = copy_col ? SVD_COPY_COL_KERNEL : SVD_COPY_ROW_KERNEL;
156  viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
158 
159  viennacl::ocl::enqueue(kernel(
160  A,
161  V,
162  static_cast<cl_uint>(row_start),
163  static_cast<cl_uint>(col_start),
164  copy_col ? static_cast<cl_uint>(A.size1())
165  : static_cast<cl_uint>(A.size2()),
166  static_cast<cl_uint>(A.internal_size2())
167  ));
168 
169  }
170 
171 
172  template<typename SCALARTYPE, unsigned int ALIGNMENT>
177  vcl_size_t row_start,
178  vcl_size_t col_start,
180  bool is_column
181  )
182  {
183  boost::numeric::ublas::vector<SCALARTYPE> tmp = boost::numeric::ublas::scalar_vector<SCALARTYPE>(size, 0);
184 
185  copy_vec(A, D, row_start, col_start, is_column);
186  fast_copy(D.begin(), D.begin() + vcl_ptrdiff_t(size - start), tmp.begin() + start);
187 
188  //std::cout << "1: " << tmp << "\n";
189 
190  detail::householder_vector(tmp, start);
191  fast_copy(tmp, D);
192 
193  //std::cout << "2: " << D << "\n";
194  }
195 
196  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename VectorType>
198  VectorType & dh,
199  VectorType & sh
200  )
201  {
204 
205  viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
207 
208  viennacl::ocl::enqueue(kernel(
209  A,
210  D,
211  S,
212  static_cast<cl_uint>(A.size1()),
213  static_cast<cl_uint>(A.size2()),
214  static_cast<cl_uint>(A.internal_size2())
215  ));
216 
217  fast_copy(D, dh);
218  fast_copy(S, sh);
219  }
220 
221  }
222  }
223 }
224 
225 #endif
viennacl::ocl::kernel & get_kernel(std::string const &program_name, std::string const &kernel_name)
Convenience function for retrieving the kernel of a program directly from the context.
Definition: context.hpp:470
size_type size1() const
Returns the number of rows.
Definition: matrix.hpp:625
std::size_t vcl_size_t
Definition: forwards.h:58
Represents an OpenCL device within ViennaCL.
const std::string SVD_COPY_ROW_KERNEL
Definition: qr-method-common.hpp:48
viennacl::ocl::kernel & get_kernel(std::string const &prog_name, std::string const &kernel_name)
Convenience function for getting the kernel for a particular program from the current active context...
Definition: backend.hpp:318
const std::string SVD_HOUSEHOLDER_UPDATE_A_RIGHT_KERNEL
Definition: qr-method-common.hpp:44
void copy_vec(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, viennacl::vector< SCALARTYPE, ALIGNMENT > &V, vcl_size_t row_start, vcl_size_t col_start, bool copy_col)
Definition: qr-method-common.hpp:147
const std::string SVD_COPY_COL_KERNEL
Definition: qr-method-common.hpp:47
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:59
Implementation of the dense matrix class.
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:51
void prepare_householder_vector(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, viennacl::vector< SCALARTYPE, ALIGNMENT > &D, vcl_size_t size, vcl_size_t row_start, vcl_size_t col_start, vcl_size_t start, bool is_column)
Definition: qr-method-common.hpp:173
void cdiv(T xr, T xi, T yr, T yi, T &cdivr, T &cdivi)
Definition: qr-method-common.hpp:124
A dense matrix class.
Definition: forwards.h:293
void normalize(VectorType &x, vcl_size_t size)
Definition: qr-method-common.hpp:87
size_type size2() const
Returns the number of columns.
Definition: matrix.hpp:627
OpenCL kernel file for singular value decomposition.
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:48
const std::string SVD_MATRIX_TRANSPOSE_KERNEL
Definition: qr-method-common.hpp:49
Main kernel class for generating OpenCL kernels for singular value decomposition of dense matrices...
Definition: svd.hpp:503
const std::string SVD_HOUSEHOLDER_UPDATE_QL_KERNEL
Definition: qr-method-common.hpp:45
const std::string SVD_FINAL_ITER_UPDATE_KERNEL
Definition: qr-method-common.hpp:53
Implementation of a smart-pointer-like class for handling OpenCL handles.
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:43
const std::string SVD_BIDIAG_PACK_KERNEL
Definition: qr-method-common.hpp:42
const std::string SVD_GIVENS_PREV_KERNEL
Definition: qr-method-common.hpp:51
void householder_vector(VectorType &v, vcl_size_t start)
Definition: qr-method-common.hpp:97
const std::string SVD_INVERSE_SIGNS_KERNEL
Definition: qr-method-common.hpp:50
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
const std::string SVD_GIVENS_NEXT_KERNEL
Definition: qr-method-common.hpp:52
void bidiag_pack(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, VectorType &dh, VectorType &sh)
Definition: qr-method-common.hpp:197
SCALARTYPE pythag(SCALARTYPE a, SCALARTYPE b)
Definition: qr-method-common.hpp:65
const std::string SVD_HOUSEHOLDER_UPDATE_A_LEFT_KERNEL
Definition: qr-method-common.hpp:43
void transpose(MatrixType &A)
Definition: qr-method-common.hpp:107
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:825
const std::string SVD_UPDATE_QR_COLUMN_KERNEL
Definition: qr-method-common.hpp:54
Representation of an OpenCL kernel in ViennaCL.
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
std::ptrdiff_t vcl_ptrdiff_t
Definition: forwards.h:59
VectorType::value_type norm_lcl(VectorType const &x, vcl_size_t size)
Definition: qr-method-common.hpp:78
void fast_copy(const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_begin, const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_end, CPU_ITERATOR cpu_begin)
STL-like transfer of a GPU vector to the CPU. The cpu type is assumed to reside in a linear piece of ...
Definition: vector.hpp:1284
A collection of compile time type deductions.
const std::string SVD_HOUSEHOLDER_UPDATE_QR_KERNEL
Definition: qr-method-common.hpp:46
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
Definition: matrix.hpp:649
SCALARTYPE sign(SCALARTYPE val)
Definition: qr-method-common.hpp:71