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;
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 = static_cast<int>(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[
vcl_size_t(l)]) <= detail::EPS)
137 goto_test_conv =
true;
141 if (std::fabs(q[
vcl_size_t(l) - 1]) <= detail::EPS)
150 CPU_ScalarType c = 0.0;
151 CPU_ScalarType s = 1.0;
156 for (
int i = l; i <= k; i++)
161 if (std::fabs(f) <= detail::EPS)
205 if (iter >= detail::ITER_MAX - 1)
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 = static_cast<vcl_size_t>(l) + 1; i <= static_cast<vcl_size_t>(k); i++)
258 givens_prev(vcl_v, tmp1, tmp2, static_cast<int>(n), l, k);
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()),
362 static_cast<cl_uint>(A.
size1()),
414 template<
typename SCALARTYPE,
unsigned int ALIGNMENT>
422 if (col_start + 1 >= A.
size2())
433 static_cast<cl_uint>(row_start),
434 static_cast<cl_uint>(col_start),
435 static_cast<cl_uint>(A.
size1()),
436 static_cast<cl_uint>(A.
size2()),
448 static_cast<cl_uint>(A.
size1()),
449 static_cast<cl_uint>(A.
size2()),
458 template<
typename SCALARTYPE,
unsigned int ALIGNMENT>
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);
523 boost::numeric::ublas::matrix<SCALARTYPE> h_Sigma(row_num, col_num);
527 h_Sigma(i, i) = dh[i];
void change_signs(MatrixType &matrix, VectorType &signs, int n)
const std::string SVD_INVERSE_SIGNS_KERNEL
Represents an OpenCL kernel within ViennaCL.
Implementation of the dense matrix class.
size_type local_work_size(int index=0) const
Returns the local work size at the respective dimension.
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
const std::string SVD_HOUSEHOLDER_UPDATE_A_RIGHT_KERNEL
const std::string SVD_HOUSEHOLDER_UPDATE_QR_KERNEL
void transpose(matrix_base< SCALARTYPE > &A)
static void init(viennacl::ocl::context &ctx)
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
T max(const T &lhs, const T &rhs)
Maximum.
void prepare_householder_vector(matrix_base< SCALARTYPE > &A, vector_base< SCALARTYPE > &D, vcl_size_t size, vcl_size_t row_start, vcl_size_t col_start, vcl_size_t start, bool is_column)
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.
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
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)
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Common routines used for the QR method and SVD. Experimental.
A class representing local (shared) OpenCL memory. Typically used as kernel argument.
Main kernel class for generating OpenCL kernels for singular value decomposition of dense matrices...
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.
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)
const std::string SVD_HOUSEHOLDER_UPDATE_A_LEFT_KERNEL
size_type size2() const
Returns the number of columns.
void bidiag_pack_svd(viennacl::matrix< SCALARTYPE > &A, VectorType &dh, VectorType &sh)
void bidiag(viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &Ai, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QL, viennacl::matrix< SCALARTYPE, row_major, ALIGNMENT > &QR)
SCALARTYPE pythag(SCALARTYPE a, SCALARTYPE b)
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
size_type size1() const
Returns the number of rows.
void givens_prev(MatrixType &matrix, VectorType &tmp1, VectorType &tmp2, int n, int l, int k)
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
const std::string SVD_HOUSEHOLDER_UPDATE_QL_KERNEL
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
size_type global_work_size(int index=0) const
Returns the global work size at the respective dimension.
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
void svd_qr_shift(MatrixType &vcl_u, MatrixType &vcl_v, CPU_VectorType &q, CPU_VectorType &e)
T min(const T &lhs, const T &rhs)
Minimum.
const std::string SVD_GIVENS_PREV_KERNEL