1 #ifndef VIENNACL_LINALG_SVD_HPP
2 #define VIENNACL_LINALG_SVD_HPP
29 #include <boost/numeric/ublas/vector.hpp>
30 #include <boost/numeric/ublas/matrix.hpp>
47 template<
typename MatrixType,
typename VectorType>
56 typedef typename MatrixType::value_type ScalarType;
69 static_cast<cl_uint>(n),
70 static_cast<cl_uint>(matrix.internal_size1()),
71 static_cast<cl_uint>(l + 1),
72 static_cast<cl_uint
>(k + 1)
77 template<
typename MatrixType,
typename VectorType>
80 typedef typename MatrixType::value_type ScalarType;
95 static_cast<cl_uint>(n),
96 static_cast<cl_uint>(matrix.internal_size1())
100 template<
typename MatrixType,
typename CPU_VectorType>
106 typedef typename MatrixType::value_type ScalarType;
109 int n =
static_cast<int>(q.size());
110 int m =
static_cast<int>(vcl_u.size1());
115 std::vector<CPU_ScalarType> signs_v(n, 1);
116 std::vector<CPU_ScalarType> cs1(n), ss1(n), cs2(n), ss2(n);
120 bool goto_test_conv =
false;
122 for (
int k = n - 1; k >= 0; k--)
127 for (iter = 0; iter < detail::ITER_MAX; iter++)
131 for (l = k; l >= 0; l--)
133 goto_test_conv =
false;
134 if (std::fabs(e[l]) <= detail::EPS)
137 goto_test_conv =
true;
141 if (std::fabs(q[l - 1]) <= detail::EPS)
150 CPU_ScalarType c = 0.0;
151 CPU_ScalarType s = 1.0;
156 for (
int i = l; i <= k; i++)
158 CPU_ScalarType f = s * e[i];
161 if (std::fabs(f) <= detail::EPS)
167 CPU_ScalarType g = q[i];
191 CPU_ScalarType z = q[k];
205 if (iter >= detail::ITER_MAX - 1)
208 CPU_ScalarType x = q[l];
209 CPU_ScalarType y = q[k - 1];
210 CPU_ScalarType g = e[k - 1];
211 CPU_ScalarType h = e[k];
212 CPU_ScalarType f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y);
214 g = detail::pythag<CPU_ScalarType>(f, 1);
217 f = ((x - z) * (x + z) + h * (y / (f - g) - h)) / x;
219 f = ((x - z) * (x + z) + h * (y / (f + g) - h)) / x;
222 CPU_ScalarType c = 1;
223 CPU_ScalarType s = 1;
225 for (
vcl_size_t i = l + 1; i <= static_cast<vcl_size_t>(k); i++)
327 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
335 if(row_start + 1 >= A.
size1())
346 static_cast<cl_uint>(row_start),
347 static_cast<cl_uint>(col_start),
348 static_cast<cl_uint>(A.
size1()),
349 static_cast<cl_uint>(A.
size2()),
361 static_cast<cl_uint>(A.
size1()),
362 static_cast<cl_uint>(A.
size2()),
412 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
420 if(col_start + 1 >= A.
size2())
431 static_cast<cl_uint>(row_start),
432 static_cast<cl_uint>(col_start),
433 static_cast<cl_uint>(A.
size1()),
434 static_cast<cl_uint>(A.
size2()),
446 static_cast<cl_uint>(A.
size1()),
447 static_cast<cl_uint>(A.
size2()),
456 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
467 vcl_size_t big_to = std::max(row_num, col_num);
491 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
514 boost::numeric::ublas::vector<SCALARTYPE> dh = boost::numeric::ublas::scalar_vector<SCALARTYPE>(to, 0);
515 boost::numeric::ublas::vector<SCALARTYPE> sh = boost::numeric::ublas::scalar_vector<SCALARTYPE>(to + 1, 0);
522 boost::numeric::ublas::matrix<SCALARTYPE> h_Sigma(row_num, col_num);
526 h_Sigma(i, i) = dh[i];
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
void change_signs(MatrixType &matrix, VectorType &signs, int n)
Definition: svd.hpp:78
const std::string SVD_HOUSEHOLDER_UPDATE_A_RIGHT_KERNEL
Definition: qr-method-common.hpp:44
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:59
Implementation of the dense matrix class.
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
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
A dense matrix class.
Definition: forwards.h:293
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
Definition: forwards.h:299
void svd(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QL, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QR)
Computes the singular value decomposition of a matrix A. Experimental in 1.3.x.
Definition: svd.hpp:492
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:245
size_type size2() const
Returns the number of columns.
Definition: matrix.hpp:627
OpenCL kernel file for singular value decomposition.
bool householder_r(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &Q, viennacl::vector< SCALARTYPE, ALIGNMENT > &D, vcl_size_t row_start, vcl_size_t col_start)
Definition: svd.hpp:413
static void init(viennacl::ocl::context &ctx)
Definition: svd.hpp:510
Common routines used for the QR method and SVD. Experimental.
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:48
A class representing local (shared) OpenCL memory. Typically used as kernel argument.
Definition: local_mem.hpp:33
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
bool householder_c(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &Q, viennacl::vector< SCALARTYPE, ALIGNMENT > &D, vcl_size_t row_start, vcl_size_t col_start)
Definition: svd.hpp:328
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
const std::string SVD_GIVENS_PREV_KERNEL
Definition: qr-method-common.hpp:51
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
void bidiag_pack(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &A, VectorType &dh, VectorType &sh)
Definition: qr-method-common.hpp:197
void bidiag(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &Ai, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QL, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QR)
Definition: svd.hpp:457
SCALARTYPE pythag(SCALARTYPE a, SCALARTYPE b)
Definition: qr-method-common.hpp:65
void givens_prev(MatrixType &matrix, VectorType &tmp1, VectorType &tmp2, int n, int l, int k)
Definition: svd.hpp:48
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
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:276
void svd_qr_shift(MatrixType &vcl_u, MatrixType &vcl_v, CPU_VectorType &q, CPU_VectorType &e)
Definition: svd.hpp:101
size_type global_work_size(int index=0) const
Returns the global work size at the respective dimension.
Definition: kernel.hpp:759
size_type local_work_size(int index=0) const
Returns the local work size at the respective dimension.
Definition: kernel.hpp:750
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
size_type size1() const
Definition: matrix.hpp:85