You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by ap...@apache.org on 2016/06/10 16:52:07 UTC
[02/51] [partial] mahout git commit: Revert "(nojira) add
native-viennaCL module to codebase. closes apache/mahout#241"
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/opencl/vector_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/opencl/vector_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/opencl/vector_operations.hpp
deleted file mode 100644
index cd04482..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/opencl/vector_operations.hpp
+++ /dev/null
@@ -1,1263 +0,0 @@
-#ifndef VIENNACL_LINALG_OPENCL_VECTOR_OPERATIONS_HPP_
-#define VIENNACL_LINALG_OPENCL_VECTOR_OPERATIONS_HPP_
-
-/* =========================================================================
- Copyright (c) 2010-2016, Institute for Microelectronics,
- Institute for Analysis and Scientific Computing,
- TU Wien.
- Portions of this software are copyright by UChicago Argonne, LLC.
-
- -----------------
- ViennaCL - The Vienna Computing Library
- -----------------
-
- Project Head: Karl Rupp rupp@iue.tuwien.ac.at
-
- (A list of authors and contributors can be found in the manual)
-
- License: MIT (X11), see file LICENSE in the base directory
-============================================================================= */
-
-/** @file viennacl/linalg/opencl/vector_operations.hpp
- @brief Implementations of vector operations using OpenCL
-*/
-
-#include <cmath>
-
-#include "viennacl/forwards.h"
-#include "viennacl/detail/vector_def.hpp"
-#include "viennacl/ocl/device.hpp"
-#include "viennacl/ocl/handle.hpp"
-#include "viennacl/ocl/kernel.hpp"
-#include "viennacl/scalar.hpp"
-#include "viennacl/tools/tools.hpp"
-#include "viennacl/linalg/opencl/common.hpp"
-#include "viennacl/linalg/opencl/kernels/vector.hpp"
-#include "viennacl/linalg/opencl/kernels/vector_element.hpp"
-#include "viennacl/linalg/opencl/kernels/scan.hpp"
-#include "viennacl/meta/predicate.hpp"
-#include "viennacl/meta/enable_if.hpp"
-#include "viennacl/traits/size.hpp"
-#include "viennacl/traits/start.hpp"
-#include "viennacl/traits/handle.hpp"
-#include "viennacl/traits/stride.hpp"
-
-namespace viennacl
-{
-namespace linalg
-{
-namespace opencl
-{
-
-//
-// Introductory note: By convention, all dimensions are already checked in the dispatcher frontend. No need to double-check again in here!
-//
-template<typename DestNumericT, typename SrcNumericT>
-void convert(vector_base<DestNumericT> & dest, vector_base<SrcNumericT> const & src)
-{
- assert(viennacl::traits::opencl_handle(dest).context() == viennacl::traits::opencl_handle(src).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- std::string kernel_name("convert_");
- kernel_name += viennacl::ocl::type_to_string<DestNumericT>::apply();
- kernel_name += "_";
- kernel_name += viennacl::ocl::type_to_string<SrcNumericT>::apply();
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(dest).context());
- viennacl::linalg::opencl::kernels::vector_convert::init(ctx);
- viennacl::ocl::kernel& k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_convert::program_name(), kernel_name);
-
- viennacl::ocl::enqueue(k( dest, cl_uint(dest.start()), cl_uint(dest.stride()), cl_uint(dest.size()),
- src, cl_uint( src.start()), cl_uint( src.stride())
- ) );
-
-}
-
-template <typename T, typename ScalarType1>
-void av(vector_base<T> & vec1,
- vector_base<T> const & vec2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- cl_uint options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(),
- (viennacl::is_cpu_scalar<ScalarType1>::value ? "av_cpu" : "av_gpu"));
- k.global_work_size(0, std::min<vcl_size_t>(128 * k.local_work_size(),
- viennacl::tools::align_to_multiple<vcl_size_t>(viennacl::traits::size(vec1), k.local_work_size()) ) );
-
- viennacl::ocl::packed_cl_uint size_vec1;
- size_vec1.start = cl_uint(viennacl::traits::start(vec1));
- size_vec1.stride = cl_uint(viennacl::traits::stride(vec1));
- size_vec1.size = cl_uint(viennacl::traits::size(vec1));
- size_vec1.internal_size = cl_uint(viennacl::traits::internal_size(vec1));
-
- viennacl::ocl::packed_cl_uint size_vec2;
- size_vec2.start = cl_uint(viennacl::traits::start(vec2));
- size_vec2.stride = cl_uint(viennacl::traits::stride(vec2));
- size_vec2.size = cl_uint(viennacl::traits::size(vec2));
- size_vec2.internal_size = cl_uint(viennacl::traits::internal_size(vec2));
-
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- size_vec1,
-
- viennacl::traits::opencl_handle(viennacl::tools::promote_if_host_scalar<T>(alpha)),
- options_alpha,
- viennacl::traits::opencl_handle(vec2),
- size_vec2 )
- );
-}
-
-
-template <typename T, typename ScalarType1, typename ScalarType2>
-void avbv(vector_base<T> & vec1,
- vector_base<T> const & vec2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
- vector_base<T> const & vec3, ScalarType2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec2).context() == viennacl::traits::opencl_handle(vec3).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- std::string kernel_name;
- if (viennacl::is_cpu_scalar<ScalarType1>::value && viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_cpu_cpu";
- else if (viennacl::is_cpu_scalar<ScalarType1>::value && !viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_cpu_gpu";
- else if (!viennacl::is_cpu_scalar<ScalarType1>::value && viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_gpu_cpu";
- else
- kernel_name = "avbv_gpu_gpu";
-
- cl_uint options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
- cl_uint options_beta = detail::make_options(len_beta, reciprocal_beta, flip_sign_beta);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), kernel_name);
- k.global_work_size(0, std::min<vcl_size_t>(128 * k.local_work_size(),
- viennacl::tools::align_to_multiple<vcl_size_t>(viennacl::traits::size(vec1), k.local_work_size()) ) );
-
- viennacl::ocl::packed_cl_uint size_vec1;
- size_vec1.start = cl_uint(viennacl::traits::start(vec1));
- size_vec1.stride = cl_uint(viennacl::traits::stride(vec1));
- size_vec1.size = cl_uint(viennacl::traits::size(vec1));
- size_vec1.internal_size = cl_uint(viennacl::traits::internal_size(vec1));
-
- viennacl::ocl::packed_cl_uint size_vec2;
- size_vec2.start = cl_uint(viennacl::traits::start(vec2));
- size_vec2.stride = cl_uint(viennacl::traits::stride(vec2));
- size_vec2.size = cl_uint(viennacl::traits::size(vec2));
- size_vec2.internal_size = cl_uint(viennacl::traits::internal_size(vec2));
-
- viennacl::ocl::packed_cl_uint size_vec3;
- size_vec3.start = cl_uint(viennacl::traits::start(vec3));
- size_vec3.stride = cl_uint(viennacl::traits::stride(vec3));
- size_vec3.size = cl_uint(viennacl::traits::size(vec3));
- size_vec3.internal_size = cl_uint(viennacl::traits::internal_size(vec3));
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- size_vec1,
-
- viennacl::traits::opencl_handle(viennacl::tools::promote_if_host_scalar<T>(alpha)),
- options_alpha,
- viennacl::traits::opencl_handle(vec2),
- size_vec2,
-
- viennacl::traits::opencl_handle(viennacl::tools::promote_if_host_scalar<T>(beta)),
- options_beta,
- viennacl::traits::opencl_handle(vec3),
- size_vec3 )
- );
-}
-
-
-template <typename T, typename ScalarType1, typename ScalarType2>
-void avbv_v(vector_base<T> & vec1,
- vector_base<T> const & vec2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
- vector_base<T> const & vec3, ScalarType2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec2).context() == viennacl::traits::opencl_handle(vec3).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- std::string kernel_name;
- if (viennacl::is_cpu_scalar<ScalarType1>::value && viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_v_cpu_cpu";
- else if (viennacl::is_cpu_scalar<ScalarType1>::value && !viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_v_cpu_gpu";
- else if (!viennacl::is_cpu_scalar<ScalarType1>::value && viennacl::is_cpu_scalar<ScalarType2>::value)
- kernel_name = "avbv_v_gpu_cpu";
- else
- kernel_name = "avbv_v_gpu_gpu";
-
- cl_uint options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
- cl_uint options_beta = detail::make_options(len_beta, reciprocal_beta, flip_sign_beta);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), kernel_name);
- k.global_work_size(0, std::min<vcl_size_t>(128 * k.local_work_size(),
- viennacl::tools::align_to_multiple<vcl_size_t>(viennacl::traits::size(vec1), k.local_work_size()) ) );
-
- viennacl::ocl::packed_cl_uint size_vec1;
- size_vec1.start = cl_uint(viennacl::traits::start(vec1));
- size_vec1.stride = cl_uint(viennacl::traits::stride(vec1));
- size_vec1.size = cl_uint(viennacl::traits::size(vec1));
- size_vec1.internal_size = cl_uint(viennacl::traits::internal_size(vec1));
-
- viennacl::ocl::packed_cl_uint size_vec2;
- size_vec2.start = cl_uint(viennacl::traits::start(vec2));
- size_vec2.stride = cl_uint(viennacl::traits::stride(vec2));
- size_vec2.size = cl_uint(viennacl::traits::size(vec2));
- size_vec2.internal_size = cl_uint(viennacl::traits::internal_size(vec2));
-
- viennacl::ocl::packed_cl_uint size_vec3;
- size_vec3.start = cl_uint(viennacl::traits::start(vec3));
- size_vec3.stride = cl_uint(viennacl::traits::stride(vec3));
- size_vec3.size = cl_uint(viennacl::traits::size(vec3));
- size_vec3.internal_size = cl_uint(viennacl::traits::internal_size(vec3));
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- size_vec1,
-
- viennacl::traits::opencl_handle(viennacl::tools::promote_if_host_scalar<T>(alpha)),
- options_alpha,
- viennacl::traits::opencl_handle(vec2),
- size_vec2,
-
- viennacl::traits::opencl_handle(viennacl::tools::promote_if_host_scalar<T>(beta)),
- options_beta,
- viennacl::traits::opencl_handle(vec3),
- size_vec3 )
- );
-}
-
-
-/** @brief Assign a constant value to a vector (-range/-slice)
-*
-* @param vec1 The vector to which the value should be assigned
-* @param alpha The value to be assigned
-* @param up_to_internal_size Specifies whether alpha should also be written to padded memory (mostly used for clearing the whole buffer).
-*/
-template <typename T>
-void vector_assign(vector_base<T> & vec1, const T & alpha, bool up_to_internal_size = false)
-{
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "assign_cpu");
- k.global_work_size(0, std::min<vcl_size_t>(128 * k.local_work_size(),
- viennacl::tools::align_to_multiple<vcl_size_t>(viennacl::traits::size(vec1), k.local_work_size()) ) );
-
- cl_uint size = up_to_internal_size ? cl_uint(vec1.internal_size()) : cl_uint(viennacl::traits::size(vec1));
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- cl_uint(viennacl::traits::start(vec1)),
- cl_uint(viennacl::traits::stride(vec1)),
- size,
- cl_uint(vec1.internal_size()), //Note: Do NOT use traits::internal_size() here, because vector proxies don't require padding.
- viennacl::traits::opencl_handle(T(alpha)) )
- );
-}
-
-
-/** @brief Swaps the contents of two vectors, data is copied
-*
-* @param vec1 The first vector (or -range, or -slice)
-* @param vec2 The second vector (or -range, or -slice)
-*/
-template <typename T>
-void vector_swap(vector_base<T> & vec1, vector_base<T> & vec2)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "swap");
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- cl_uint(viennacl::traits::start(vec1)),
- cl_uint(viennacl::traits::stride(vec1)),
- cl_uint(viennacl::traits::size(vec1)),
- viennacl::traits::opencl_handle(vec2),
- cl_uint(viennacl::traits::start(vec2)),
- cl_uint(viennacl::traits::stride(vec2)),
- cl_uint(viennacl::traits::size(vec2)))
- );
-}
-
-///////////////////////// Binary Elementwise operations /////////////
-
-/** @brief Implementation of the element-wise operation v1 = v2 .* v3 and v1 = v2 ./ v3 (using MATLAB syntax)
-*
-* @param vec1 The result vector (or -range, or -slice)
-* @param proxy The proxy object holding v2, v3 and the operation
-*/
-template <typename T, typename OP>
-void element_op(vector_base<T> & vec1,
- vector_expression<const vector_base<T>, const vector_base<T>, op_element_binary<OP> > const & proxy)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(proxy.lhs()).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(proxy.rhs()).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector_element<T>::init(ctx);
-
- std::string kernel_name = "element_pow";
- cl_uint op_type = 2; //0: product, 1: division, 2: power
- if (viennacl::is_division<OP>::value)
- {
- op_type = 1;
- kernel_name = "element_div";
- }
- else if (viennacl::is_product<OP>::value)
- {
- op_type = 0;
- kernel_name = "element_prod";
- }
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_element<T>::program_name(), kernel_name);
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- cl_uint(viennacl::traits::start(vec1)),
- cl_uint(viennacl::traits::stride(vec1)),
- cl_uint(viennacl::traits::size(vec1)),
-
- viennacl::traits::opencl_handle(proxy.lhs()),
- cl_uint(viennacl::traits::start(proxy.lhs())),
- cl_uint(viennacl::traits::stride(proxy.lhs())),
-
- viennacl::traits::opencl_handle(proxy.rhs()),
- cl_uint(viennacl::traits::start(proxy.rhs())),
- cl_uint(viennacl::traits::stride(proxy.rhs())),
-
- op_type)
- );
-}
-
-///////////////////////// Unary Elementwise operations /////////////
-
-/** @brief Implementation of unary element-wise operations v1 = OP(v2)
-*
-* @param vec1 The result vector (or -range, or -slice)
-* @param proxy The proxy object holding v2 and the operation
-*/
-template <typename T, typename OP>
-void element_op(vector_base<T> & vec1,
- vector_expression<const vector_base<T>, const vector_base<T>, op_element_unary<OP> > const & proxy)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(proxy.lhs()).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(proxy.rhs()).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector_element<T>::init(ctx);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_element<T>::program_name(), detail::op_to_string(OP()) + "_assign");
-
- viennacl::ocl::packed_cl_uint size_vec1;
- size_vec1.start = cl_uint(viennacl::traits::start(vec1));
- size_vec1.stride = cl_uint(viennacl::traits::stride(vec1));
- size_vec1.size = cl_uint(viennacl::traits::size(vec1));
- size_vec1.internal_size = cl_uint(viennacl::traits::internal_size(vec1));
-
- viennacl::ocl::packed_cl_uint size_vec2;
- size_vec2.start = cl_uint(viennacl::traits::start(proxy.lhs()));
- size_vec2.stride = cl_uint(viennacl::traits::stride(proxy.lhs()));
- size_vec2.size = cl_uint(viennacl::traits::size(proxy.lhs()));
- size_vec2.internal_size = cl_uint(viennacl::traits::internal_size(proxy.lhs()));
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- size_vec1,
- viennacl::traits::opencl_handle(proxy.lhs()),
- size_vec2)
- );
-}
-
-///////////////////////// Norms and inner product ///////////////////
-
-/** @brief Computes the partial inner product of two vectors - implementation. Library users should call inner_prod(vec1, vec2).
-*
-* @param vec1 The first vector
-* @param vec2 The second vector
-* @param partial_result The results of each group
-*/
-template <typename T>
-void inner_prod_impl(vector_base<T> const & vec1,
- vector_base<T> const & vec2,
- vector_base<T> & partial_result)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec2).context() == viennacl::traits::opencl_handle(partial_result).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- assert( (viennacl::traits::size(vec1) == viennacl::traits::size(vec2))
- && bool("Incompatible vector sizes in inner_prod_impl()!"));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "inner_prod1");
-
- assert( (k.global_work_size() / k.local_work_size() <= partial_result.size()) && bool("Size mismatch for partial reduction in inner_prod_impl()") );
-
- viennacl::ocl::packed_cl_uint size_vec1;
- size_vec1.start = cl_uint(viennacl::traits::start(vec1));
- size_vec1.stride = cl_uint(viennacl::traits::stride(vec1));
- size_vec1.size = cl_uint(viennacl::traits::size(vec1));
- size_vec1.internal_size = cl_uint(viennacl::traits::internal_size(vec1));
-
- viennacl::ocl::packed_cl_uint size_vec2;
- size_vec2.start = cl_uint(viennacl::traits::start(vec2));
- size_vec2.stride = cl_uint(viennacl::traits::stride(vec2));
- size_vec2.size = cl_uint(viennacl::traits::size(vec2));
- size_vec2.internal_size = cl_uint(viennacl::traits::internal_size(vec2));
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- size_vec1,
- viennacl::traits::opencl_handle(vec2),
- size_vec2,
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(partial_result)
- )
- );
-}
-
-
-//implementation of inner product:
-//namespace {
-/** @brief Computes the inner product of two vectors - implementation. Library users should call inner_prod(vec1, vec2).
-*
-* @param vec1 The first vector
-* @param vec2 The second vector
-* @param result The result scalar (on the gpu)
-*/
-template <typename T>
-void inner_prod_impl(vector_base<T> const & vec1,
- vector_base<T> const & vec2,
- scalar<T> & result)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
-
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec1));
- temp.resize(work_groups, ctx); // bring default-constructed vectors to the correct size:
-
- // Step 1: Compute partial inner products for each work group:
- inner_prod_impl(vec1, vec2, temp);
-
- // Step 2: Sum partial results:
- viennacl::ocl::kernel & ksum = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "sum");
-
- ksum.global_work_size(0, ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- cl_uint(1),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result) )
- );
-}
-
-namespace detail
-{
- template<typename NumericT>
- viennacl::ocl::packed_cl_uint make_layout(vector_base<NumericT> const & vec)
- {
- viennacl::ocl::packed_cl_uint ret;
- ret.start = cl_uint(viennacl::traits::start(vec));
- ret.stride = cl_uint(viennacl::traits::stride(vec));
- ret.size = cl_uint(viennacl::traits::size(vec));
- ret.internal_size = cl_uint(viennacl::traits::internal_size(vec));
- return ret;
- }
-}
-
-/** @brief Computes multiple inner products where one argument is common to all inner products. <x, y1>, <x, y2>, ..., <x, yN>
-*
-* @param x The common vector
-* @param vec_tuple The tuple of vectors y1, y2, ..., yN
-* @param result The result vector
-*/
-template <typename NumericT>
-void inner_prod_impl(vector_base<NumericT> const & x,
- vector_tuple<NumericT> const & vec_tuple,
- vector_base<NumericT> & result)
-{
- assert(viennacl::traits::opencl_handle(x).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
- viennacl::linalg::opencl::kernels::vector<NumericT>::init(ctx);
- viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::init(ctx);
-
- viennacl::ocl::packed_cl_uint layout_x = detail::make_layout(x);
-
- viennacl::ocl::kernel & ksum = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::program_name(), "sum_inner_prod");
- viennacl::ocl::kernel & inner_prod_kernel_1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<NumericT>::program_name(), "inner_prod1");
- viennacl::ocl::kernel & inner_prod_kernel_2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::program_name(), "inner_prod2");
- viennacl::ocl::kernel & inner_prod_kernel_3 = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::program_name(), "inner_prod3");
- viennacl::ocl::kernel & inner_prod_kernel_4 = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::program_name(), "inner_prod4");
- viennacl::ocl::kernel & inner_prod_kernel_8 = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector_multi_inner_prod<NumericT>::program_name(), "inner_prod8");
-
- vcl_size_t work_groups = inner_prod_kernel_8.global_work_size(0) / inner_prod_kernel_8.local_work_size(0);
- viennacl::vector<NumericT> temp(8 * work_groups, viennacl::traits::context(x));
-
- vcl_size_t current_index = 0;
- while (current_index < vec_tuple.const_size())
- {
- switch (vec_tuple.const_size() - current_index)
- {
- case 7:
- case 6:
- case 5:
- case 4:
- {
- vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index );
- vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
- vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
- vector_base<NumericT> const & y3 = vec_tuple.const_at(current_index + 3);
- viennacl::ocl::enqueue(inner_prod_kernel_4( viennacl::traits::opencl_handle(x), layout_x,
- viennacl::traits::opencl_handle(y0), detail::make_layout(y0),
- viennacl::traits::opencl_handle(y1), detail::make_layout(y1),
- viennacl::traits::opencl_handle(y2), detail::make_layout(y2),
- viennacl::traits::opencl_handle(y3), detail::make_layout(y3),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 4 * inner_prod_kernel_4.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ) );
-
- ksum.global_work_size(0, 4 * ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(work_groups),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 4 * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result),
- cl_uint(viennacl::traits::start(result) + current_index * viennacl::traits::stride(result)),
- cl_uint(viennacl::traits::stride(result))
- )
- );
- }
- current_index += 4;
- break;
-
- case 3:
- {
- vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index );
- vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
- vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
- viennacl::ocl::enqueue(inner_prod_kernel_3( viennacl::traits::opencl_handle(x), layout_x,
- viennacl::traits::opencl_handle(y0), detail::make_layout(y0),
- viennacl::traits::opencl_handle(y1), detail::make_layout(y1),
- viennacl::traits::opencl_handle(y2), detail::make_layout(y2),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 3 * inner_prod_kernel_3.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ) );
-
- ksum.global_work_size(0, 3 * ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(work_groups),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 3 * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result),
- cl_uint(viennacl::traits::start(result) + current_index * viennacl::traits::stride(result)),
- cl_uint(viennacl::traits::stride(result))
- )
- );
- }
- current_index += 3;
- break;
-
- case 2:
- {
- vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index );
- vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
- viennacl::ocl::enqueue(inner_prod_kernel_2( viennacl::traits::opencl_handle(x), layout_x,
- viennacl::traits::opencl_handle(y0), detail::make_layout(y0),
- viennacl::traits::opencl_handle(y1), detail::make_layout(y1),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 2 * inner_prod_kernel_2.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ) );
-
- ksum.global_work_size(0, 2 * ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(work_groups),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 2 * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result),
- cl_uint(viennacl::traits::start(result) + current_index * viennacl::traits::stride(result)),
- cl_uint(viennacl::traits::stride(result))
- )
- );
- }
- current_index += 2;
- break;
-
- case 1:
- {
- vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index );
- viennacl::ocl::enqueue(inner_prod_kernel_1( viennacl::traits::opencl_handle(x), layout_x,
- viennacl::traits::opencl_handle(y0), detail::make_layout(y0),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 1 * inner_prod_kernel_1.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ) );
-
- ksum.global_work_size(0, 1 * ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(work_groups),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 1 * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result),
- cl_uint(viennacl::traits::start(result) + current_index * viennacl::traits::stride(result)),
- cl_uint(viennacl::traits::stride(result))
- )
- );
- }
- current_index += 1;
- break;
-
- default: //8 or more vectors
- {
- vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index );
- vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
- vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
- vector_base<NumericT> const & y3 = vec_tuple.const_at(current_index + 3);
- vector_base<NumericT> const & y4 = vec_tuple.const_at(current_index + 4);
- vector_base<NumericT> const & y5 = vec_tuple.const_at(current_index + 5);
- vector_base<NumericT> const & y6 = vec_tuple.const_at(current_index + 6);
- vector_base<NumericT> const & y7 = vec_tuple.const_at(current_index + 7);
- viennacl::ocl::enqueue(inner_prod_kernel_8( viennacl::traits::opencl_handle(x), layout_x,
- viennacl::traits::opencl_handle(y0), detail::make_layout(y0),
- viennacl::traits::opencl_handle(y1), detail::make_layout(y1),
- viennacl::traits::opencl_handle(y2), detail::make_layout(y2),
- viennacl::traits::opencl_handle(y3), detail::make_layout(y3),
- viennacl::traits::opencl_handle(y4), detail::make_layout(y4),
- viennacl::traits::opencl_handle(y5), detail::make_layout(y5),
- viennacl::traits::opencl_handle(y6), detail::make_layout(y6),
- viennacl::traits::opencl_handle(y7), detail::make_layout(y7),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 8 * inner_prod_kernel_8.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ) );
-
- ksum.global_work_size(0, 8 * ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(work_groups),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * 8 * ksum.local_work_size()),
- viennacl::traits::opencl_handle(result),
- cl_uint(viennacl::traits::start(result) + current_index * viennacl::traits::stride(result)),
- cl_uint(viennacl::traits::stride(result))
- )
- );
- }
- current_index += 8;
- break;
- }
- }
-
-}
-
-
-
-//implementation of inner product:
-//namespace {
-/** @brief Computes the inner product of two vectors - implementation. Library users should call inner_prod(vec1, vec2).
-*
-* @param vec1 The first vector
-* @param vec2 The second vector
-* @param result The result scalar (on the gpu)
-*/
-template <typename T>
-void inner_prod_cpu(vector_base<T> const & vec1,
- vector_base<T> const & vec2,
- T & result)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Vectors do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
-
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec1));
- temp.resize(work_groups, ctx); // bring default-constructed vectors to the correct size:
-
- // Step 1: Compute partial inner products for each work group:
- inner_prod_impl(vec1, vec2, temp);
-
- // Step 2: Sum partial results:
-
- // Now copy partial results from GPU back to CPU and run reduction there:
- std::vector<T> temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = 0;
- for (typename std::vector<T>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result += *it;
-}
-
-
-//////////// Helper for norms
-
-/** @brief Computes the partial work group results for vector norms
-*
-* @param vec The vector
-* @param partial_result The result scalar
-* @param norm_id Norm selector. 0: norm_inf, 1: norm_1, 2: norm_2
-*/
-template <typename T>
-void norm_reduction_impl(vector_base<T> const & vec,
- vector_base<T> & partial_result,
- cl_uint norm_id)
-{
- assert(viennacl::traits::opencl_handle(vec).context() == viennacl::traits::opencl_handle(partial_result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "norm");
-
- assert( (k.global_work_size() / k.local_work_size() <= partial_result.size()) && bool("Size mismatch for partial reduction in norm_reduction_impl()") );
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec),
- cl_uint(viennacl::traits::start(vec)),
- cl_uint(viennacl::traits::stride(vec)),
- cl_uint(viennacl::traits::size(vec)),
- cl_uint(norm_id),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(partial_result) )
- );
-}
-
-
-//////////// Norm 1
-
-/** @brief Computes the l^1-norm of a vector
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_1_impl(vector_base<T> const & vec,
- scalar<T> & result)
-{
- assert(viennacl::traits::opencl_handle(vec).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec).context());
-
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 1);
-
- // Step 2: Compute the partial reduction using OpenCL
- viennacl::ocl::kernel & ksum = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "sum");
-
- ksum.global_work_size(0, ksum.local_work_size(0));
- viennacl::ocl::enqueue(ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- cl_uint(1),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * ksum.local_work_size()),
- result)
- );
-}
-
-/** @brief Computes the l^1-norm of a vector with final reduction on CPU
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_1_cpu(vector_base<T> const & vec,
- T & result)
-{
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 1);
-
- // Step 2: Now copy partial results from GPU back to CPU and run reduction there:
- typedef std::vector<typename viennacl::result_of::cl_type<T>::type> CPUVectorType;
-
- CPUVectorType temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = 0;
- for (typename CPUVectorType::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result += static_cast<T>(*it);
-}
-
-
-
-//////// Norm 2
-
-
-/** @brief Computes the l^2-norm of a vector - implementation using OpenCL summation at second step
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_2_impl(vector_base<T> const & vec,
- scalar<T> & result)
-{
- assert(viennacl::traits::opencl_handle(vec).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec).context());
-
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 2);
-
- // Step 2: Reduction via OpenCL
- viennacl::ocl::kernel & ksum = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "sum");
-
- ksum.global_work_size(0, ksum.local_work_size(0));
- viennacl::ocl::enqueue( ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- cl_uint(2),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * ksum.local_work_size()),
- result)
- );
-}
-
-/** @brief Computes the l^1-norm of a vector with final reduction on CPU
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_2_cpu(vector_base<T> const & vec,
- T & result)
-{
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 2);
-
- // Step 2: Now copy partial results from GPU back to CPU and run reduction there:
- typedef std::vector<typename viennacl::result_of::cl_type<T>::type> CPUVectorType;
-
- CPUVectorType temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = 0;
- for (typename CPUVectorType::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result += static_cast<T>(*it);
- result = std::sqrt(result);
-}
-
-
-
-////////// Norm inf
-
-/** @brief Computes the supremum-norm of a vector
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_inf_impl(vector_base<T> const & vec,
- scalar<T> & result)
-{
- assert(viennacl::traits::opencl_handle(vec).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec).context());
-
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 0);
-
- //part 2: parallel reduction of reduced kernel:
- viennacl::ocl::kernel & ksum = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "sum");
-
- ksum.global_work_size(0, ksum.local_work_size(0));
- viennacl::ocl::enqueue( ksum(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- cl_uint(0),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * ksum.local_work_size()),
- result)
- );
-}
-
-/** @brief Computes the supremum-norm of a vector
-*
-* @param vec The vector
-* @param result The result scalar
-*/
-template <typename T>
-void norm_inf_cpu(vector_base<T> const & vec,
- T & result)
-{
- vcl_size_t work_groups = 128;
- viennacl::vector<T> temp(work_groups, viennacl::traits::context(vec));
-
- // Step 1: Compute the partial work group results
- norm_reduction_impl(vec, temp, 0);
-
- // Step 2: Now copy partial results from GPU back to CPU and run reduction there:
- typedef std::vector<typename viennacl::result_of::cl_type<T>::type> CPUVectorType;
-
- CPUVectorType temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = 0;
- for (typename CPUVectorType::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result = std::max(result, static_cast<T>(*it));
-}
-
-
-/////////// index norm_inf
-
-//This function should return a CPU scalar, otherwise statements like
-// vcl_rhs[index_norm_inf(vcl_rhs)]
-// are ambiguous
-/** @brief Computes the index of the first entry that is equal to the supremum-norm in modulus.
-*
-* @param vec The vector
-* @return The result. Note that the result must be a CPU scalar (unsigned int), since gpu scalars are floating point types.
-*/
-template <typename T>
-cl_uint index_norm_inf(vector_base<T> const & vec)
-{
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- viennacl::ocl::handle<cl_mem> h = ctx.create_memory(CL_MEM_READ_WRITE, sizeof(cl_uint));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "index_norm_inf");
- //cl_uint size = static_cast<cl_uint>(vcl_vec.internal_size());
-
- //TODO: Use multi-group kernel for large vector sizes
-
- k.global_work_size(0, k.local_work_size());
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec),
- cl_uint(viennacl::traits::start(vec)),
- cl_uint(viennacl::traits::stride(vec)),
- cl_uint(viennacl::traits::size(vec)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<T>::type) * k.local_work_size()),
- viennacl::ocl::local_mem(sizeof(cl_uint) * k.local_work_size()), h));
-
- //read value:
- cl_uint result;
- cl_int err = clEnqueueReadBuffer(ctx.get_queue().handle().get(), h.get(), CL_TRUE, 0, sizeof(cl_uint), &result, 0, NULL, NULL);
- VIENNACL_ERR_CHECK(err);
- return result;
-}
-
-
-////////// max
-
-/** @brief Computes the maximum value of a vector, where the result is stored in an OpenCL buffer.
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void max_impl(vector_base<NumericT> const & x,
- scalar<NumericT> & result)
-{
- assert(viennacl::traits::opencl_handle(x).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
- viennacl::linalg::opencl::kernels::vector<NumericT>::init(ctx);
-
- vcl_size_t work_groups = 128;
- viennacl::vector<NumericT> temp(work_groups, viennacl::traits::context(x));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<NumericT>::program_name(), "max_kernel");
-
- k.global_work_size(0, work_groups * k.local_work_size(0));
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(x),
- cl_uint(viennacl::traits::start(x)),
- cl_uint(viennacl::traits::stride(x)),
- cl_uint(viennacl::traits::size(x)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ));
-
- k.global_work_size(0, k.local_work_size());
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(result)
- ));
-}
-
-/** @brief Computes the maximum value of a vector, where the value is stored in a host value.
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void max_cpu(vector_base<NumericT> const & x,
- NumericT & result)
-{
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
- viennacl::linalg::opencl::kernels::vector<NumericT>::init(ctx);
-
- vcl_size_t work_groups = 128;
- viennacl::vector<NumericT> temp(work_groups, viennacl::traits::context(x));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<NumericT>::program_name(), "max_kernel");
-
- k.global_work_size(0, work_groups * k.local_work_size(0));
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(x),
- cl_uint(viennacl::traits::start(x)),
- cl_uint(viennacl::traits::stride(x)),
- cl_uint(viennacl::traits::size(x)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ));
-
- // Step 2: Now copy partial results from GPU back to CPU and run reduction there:
- typedef std::vector<typename viennacl::result_of::cl_type<NumericT>::type> CPUVectorType;
-
- CPUVectorType temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = static_cast<NumericT>(temp_cpu[0]);
- for (typename CPUVectorType::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result = std::max(result, static_cast<NumericT>(*it));
-
-}
-
-
-////////// min
-
-/** @brief Computes the minimum of a vector, where the result is stored in an OpenCL buffer.
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void min_impl(vector_base<NumericT> const & x,
- scalar<NumericT> & result)
-{
- assert(viennacl::traits::opencl_handle(x).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
- viennacl::linalg::opencl::kernels::vector<NumericT>::init(ctx);
-
- vcl_size_t work_groups = 128;
- viennacl::vector<NumericT> temp(work_groups, viennacl::traits::context(x));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<NumericT>::program_name(), "min_kernel");
-
- k.global_work_size(0, work_groups * k.local_work_size(0));
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(x),
- cl_uint(viennacl::traits::start(x)),
- cl_uint(viennacl::traits::stride(x)),
- cl_uint(viennacl::traits::size(x)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ));
-
- k.global_work_size(0, k.local_work_size());
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(temp),
- cl_uint(viennacl::traits::start(temp)),
- cl_uint(viennacl::traits::stride(temp)),
- cl_uint(viennacl::traits::size(temp)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(result)
- ));
-}
-
-/** @brief Computes the minimum of a vector, where the result is stored on a CPU scalar.
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void min_cpu(vector_base<NumericT> const & x,
- NumericT & result)
-{
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
- viennacl::linalg::opencl::kernels::vector<NumericT>::init(ctx);
-
- vcl_size_t work_groups = 128;
- viennacl::vector<NumericT> temp(work_groups, viennacl::traits::context(x));
-
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<NumericT>::program_name(), "min_kernel");
-
- k.global_work_size(0, work_groups * k.local_work_size(0));
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(x),
- cl_uint(viennacl::traits::start(x)),
- cl_uint(viennacl::traits::stride(x)),
- cl_uint(viennacl::traits::size(x)),
- viennacl::ocl::local_mem(sizeof(typename viennacl::result_of::cl_type<NumericT>::type) * k.local_work_size()),
- viennacl::traits::opencl_handle(temp)
- ));
-
- // Step 2: Now copy partial results from GPU back to CPU and run reduction there:
- typedef std::vector<typename viennacl::result_of::cl_type<NumericT>::type> CPUVectorType;
-
- CPUVectorType temp_cpu(work_groups);
- viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
-
- result = static_cast<NumericT>(temp_cpu[0]);
- for (typename CPUVectorType::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
- result = std::min(result, static_cast<NumericT>(*it));
-}
-
-////////// sum
-
-/** @brief Computes the sum over all entries of a vector
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void sum_impl(vector_base<NumericT> const & x,
- scalar<NumericT> & result)
-{
- assert(viennacl::traits::opencl_handle(x).context() == viennacl::traits::opencl_handle(result).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::vector<NumericT> all_ones = viennacl::scalar_vector<NumericT>(x.size(), NumericT(1), viennacl::traits::context(x));
- viennacl::linalg::opencl::inner_prod_impl(x, all_ones, result);
-}
-
-/** @brief Computes the sum over all entries of a vector.
-*
-* @param x The vector
-* @param result The result scalar
-*/
-template<typename NumericT>
-void sum_cpu(vector_base<NumericT> const & x, NumericT & result)
-{
- scalar<NumericT> tmp(0, viennacl::traits::context(x));
- sum_impl(x, tmp);
- result = tmp;
-}
-
-
-//TODO: Special case vec1 == vec2 allows improvement!!
-/** @brief Computes a plane rotation of two vectors.
-*
-* Computes (x,y) <- (alpha * x + beta * y, -beta * x + alpha * y)
-*
-* @param vec1 The first vector
-* @param vec2 The second vector
-* @param alpha The first transformation coefficient
-* @param beta The second transformation coefficient
-*/
-template <typename T>
-void plane_rotation(vector_base<T> & vec1,
- vector_base<T> & vec2,
- T alpha, T beta)
-{
- assert(viennacl::traits::opencl_handle(vec1).context() == viennacl::traits::opencl_handle(vec2).context() && bool("Operands do not reside in the same OpenCL context. Automatic migration not yet supported!"));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(vec1).context());
- viennacl::linalg::opencl::kernels::vector<T>::init(ctx);
-
- assert(viennacl::traits::size(vec1) == viennacl::traits::size(vec2));
- viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::vector<T>::program_name(), "plane_rotation");
-
- viennacl::ocl::enqueue(k(viennacl::traits::opencl_handle(vec1),
- cl_uint(viennacl::traits::start(vec1)),
- cl_uint(viennacl::traits::stride(vec1)),
- cl_uint(viennacl::traits::size(vec1)),
- viennacl::traits::opencl_handle(vec2),
- cl_uint(viennacl::traits::start(vec2)),
- cl_uint(viennacl::traits::stride(vec2)),
- cl_uint(viennacl::traits::size(vec2)),
- viennacl::traits::opencl_handle(alpha),
- viennacl::traits::opencl_handle(beta))
- );
-}
-
-
-//////////////////////////
-
-
-namespace detail
-{
- /** @brief Worker routine for scan routines using OpenCL
- *
- * Note on performance: For non-in-place scans one could optimize away the temporary 'opencl_carries'-array.
- * This, however, only provides small savings in the latency-dominated regime, yet would effectively double the amount of code to maintain.
- */
- template<typename NumericT>
- void scan_impl(vector_base<NumericT> const & input,
- vector_base<NumericT> & output,
- bool is_inclusive)
- {
- vcl_size_t local_worksize = 128;
- vcl_size_t workgroups = 128;
-
- viennacl::backend::mem_handle opencl_carries;
- viennacl::backend::memory_create(opencl_carries, sizeof(NumericT)*workgroups, viennacl::traits::context(input));
-
- viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(input).context());
- viennacl::linalg::opencl::kernels::scan<NumericT>::init(ctx);
- viennacl::ocl::kernel& k1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::scan<NumericT>::program_name(), "scan_1");
- viennacl::ocl::kernel& k2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::scan<NumericT>::program_name(), "scan_2");
- viennacl::ocl::kernel& k3 = ctx.get_kernel(viennacl::linalg::opencl::kernels::scan<NumericT>::program_name(), "scan_3");
-
- // First step: Scan within each thread group and write carries
- k1.local_work_size(0, local_worksize);
- k1.global_work_size(0, workgroups * local_worksize);
- viennacl::ocl::enqueue(k1( input, cl_uint( input.start()), cl_uint( input.stride()), cl_uint(input.size()),
- output, cl_uint(output.start()), cl_uint(output.stride()),
- cl_uint(is_inclusive ? 0 : 1), opencl_carries.opencl_handle())
- );
-
- // Second step: Compute offset for each thread group (exclusive scan for each thread group)
- k2.local_work_size(0, workgroups);
- k2.global_work_size(0, workgroups);
- viennacl::ocl::enqueue(k2(opencl_carries.opencl_handle()));
-
- // Third step: Offset each thread group accordingly
- k3.local_work_size(0, local_worksize);
- k3.global_work_size(0, workgroups * local_worksize);
- viennacl::ocl::enqueue(k3(output, cl_uint(output.start()), cl_uint(output.stride()), cl_uint(output.size()),
- opencl_carries.opencl_handle())
- );
- }
-}
-
-
-/** @brief This function implements an inclusive scan using CUDA.
-*
-* @param input Input vector.
-* @param output The output vector. Either idential to input or non-overlapping.
-*/
-template<typename NumericT>
-void inclusive_scan(vector_base<NumericT> const & input,
- vector_base<NumericT> & output)
-{
- detail::scan_impl(input, output, true);
-}
-
-
-/** @brief This function implements an exclusive scan using CUDA.
-*
-* @param input Input vector
-* @param output The output vector. Either idential to input or non-overlapping.
-*/
-template<typename NumericT>
-void exclusive_scan(vector_base<NumericT> const & input,
- vector_base<NumericT> & output)
-{
- detail::scan_impl(input, output, false);
-}
-
-
-} //namespace opencl
-} //namespace linalg
-} //namespace viennacl
-
-
-#endif
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/power_iter.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/power_iter.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/power_iter.hpp
deleted file mode 100644
index 9721517..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/power_iter.hpp
+++ /dev/null
@@ -1,129 +0,0 @@
-#ifndef VIENNACL_LINALG_POWER_ITER_HPP_
-#define VIENNACL_LINALG_POWER_ITER_HPP_
-
-/* =========================================================================
- Copyright (c) 2010-2016, Institute for Microelectronics,
- Institute for Analysis and Scientific Computing,
- TU Wien.
- Portions of this software are copyright by UChicago Argonne, LLC.
-
- -----------------
- ViennaCL - The Vienna Computing Library
- -----------------
-
- Project Head: Karl Rupp rupp@iue.tuwien.ac.at
-
- (A list of authors and contributors can be found in the manual)
-
- License: MIT (X11), see file LICENSE in the base directory
-============================================================================= */
-
-/** @file viennacl/linalg/power_iter.hpp
- @brief Defines a tag for the configuration of the power iteration method.
-
- Contributed by Astrid Rupp.
-*/
-
-#include <cmath>
-#include <vector>
-#include "viennacl/linalg/bisect.hpp"
-#include "viennacl/linalg/prod.hpp"
-#include "viennacl/linalg/norm_2.hpp"
-
-namespace viennacl
-{
- namespace linalg
- {
- /** @brief A tag for the power iteration algorithm. */
- class power_iter_tag
- {
- public:
-
- /** @brief The constructor
- *
- * @param tfac If the eigenvalue does not change more than this termination factor, the algorithm stops
- * @param max_iters Maximum number of iterations for the power iteration
- */
- power_iter_tag(double tfac = 1e-8, vcl_size_t max_iters = 50000) : termination_factor_(tfac), max_iterations_(max_iters) {}
-
- /** @brief Sets the factor for termination */
- void factor(double fct){ termination_factor_ = fct; }
-
- /** @brief Returns the factor for termination */
- double factor() const { return termination_factor_; }
-
- vcl_size_t max_iterations() const { return max_iterations_; }
- void max_iterations(vcl_size_t new_max) { max_iterations_ = new_max; }
-
- private:
- double termination_factor_;
- vcl_size_t max_iterations_;
-
- };
-
- /**
- * @brief Implementation of the calculation of the largest eigenvalue (in modulus) and the associated eigenvector using power iteration
- *
- * @param A The system matrix
- * @param tag Tag with termination factor
- * @param eigenvec Vector which holds the associated eigenvector once the routine completes
- * @return Returns the largest eigenvalue computed by the power iteration method
- */
- template<typename MatrixT, typename VectorT >
- typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
- eig(MatrixT const& A, power_iter_tag const & tag, VectorT & eigenvec)
- {
-
- typedef typename viennacl::result_of::value_type<MatrixT>::type ScalarType;
- typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type CPU_ScalarType;
-
- vcl_size_t matrix_size = A.size1();
- VectorT r(eigenvec);
- std::vector<CPU_ScalarType> s(matrix_size);
-
- for (vcl_size_t i=0; i<s.size(); ++i)
- s[i] = CPU_ScalarType(i % 3) * CPU_ScalarType(0.1234) - CPU_ScalarType(0.5); //'random' starting vector
-
- detail::copy_vec_to_vec(s, eigenvec);
-
- double epsilon = tag.factor();
- CPU_ScalarType norm = norm_2(eigenvec);
- CPU_ScalarType norm_prev = 0;
- long numiter = 0;
-
- for (vcl_size_t i=0; i<tag.max_iterations(); ++i)
- {
- if (std::fabs(norm - norm_prev) / std::fabs(norm) < epsilon)
- break;
-
- eigenvec /= norm;
- r = viennacl::linalg::prod(A, eigenvec); //using helper vector r for the computation of x <- A * x in order to avoid the repeated creation of temporaries
- eigenvec = r;
- norm_prev = norm;
- norm = norm_2(eigenvec);
- numiter++;
- }
-
- return norm;
- }
-
- /**
- * @brief Implementation of the calculation of eigenvalues using power iteration. Does not return the eigenvector.
- *
- * @param A The system matrix
- * @param tag Tag with termination factor
- * @return Returns the largest eigenvalue computed by the power iteration method
- */
- template< typename MatrixT >
- typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
- eig(MatrixT const& A, power_iter_tag const & tag)
- {
- typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type VectorT;
-
- VectorT eigenvec(A.size1());
- return eig(A, tag, eigenvec);
- }
-
- } // end namespace linalg
-} // end namespace viennacl
-#endif
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/prod.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/prod.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/prod.hpp
deleted file mode 100644
index af041dc..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/prod.hpp
+++ /dev/null
@@ -1,370 +0,0 @@
-#ifndef VIENNACL_LINALG_PROD_HPP_
-#define VIENNACL_LINALG_PROD_HPP_
-
-/* =========================================================================
- Copyright (c) 2010-2016, Institute for Microelectronics,
- Institute for Analysis and Scientific Computing,
- TU Wien.
- Portions of this software are copyright by UChicago Argonne, LLC.
-
- -----------------
- ViennaCL - The Vienna Computing Library
- -----------------
-
- Project Head: Karl Rupp rupp@iue.tuwien.ac.at
-
- (A list of authors and contributors can be found in the manual)
-
- License: MIT (X11), see file LICENSE in the base directory
-============================================================================= */
-
-/** @file viennacl/linalg/prod.hpp
- @brief Generic interface for matrix-vector and matrix-matrix products.
- See viennacl/linalg/vector_operations.hpp, viennacl/linalg/matrix_operations.hpp, and
- viennacl/linalg/sparse_matrix_operations.hpp for implementations.
-*/
-
-#include "viennacl/forwards.h"
-#include "viennacl/tools/tools.hpp"
-#include "viennacl/meta/enable_if.hpp"
-#include "viennacl/meta/tag_of.hpp"
-#include <vector>
-#include <map>
-
-namespace viennacl
-{
- //
- // generic prod function
- // uses tag dispatch to identify which algorithm
- // should be called
- //
- namespace linalg
- {
- #ifdef VIENNACL_WITH_MTL4
- // ----------------------------------------------------
- // mtl4
- //
- template< typename MatrixT, typename VectorT >
- typename viennacl::enable_if< viennacl::is_mtl4< typename viennacl::traits::tag_of< MatrixT >::type >::value,
- VectorT>::type
- prod(MatrixT const& matrix, VectorT const& vector)
- {
- return VectorT(matrix * vector);
- }
- #endif
-
- #ifdef VIENNACL_WITH_ARMADILLO
- // ----------------------------------------------------
- // Armadillo
- //
- template<typename NumericT, typename VectorT>
- VectorT prod(arma::SpMat<NumericT> const& A, VectorT const& vector)
- {
- return A * vector;
- }
- #endif
-
- #ifdef VIENNACL_WITH_EIGEN
- // ----------------------------------------------------
- // Eigen
- //
- template< typename MatrixT, typename VectorT >
- typename viennacl::enable_if< viennacl::is_eigen< typename viennacl::traits::tag_of< MatrixT >::type >::value,
- VectorT>::type
- prod(MatrixT const& matrix, VectorT const& vector)
- {
- return matrix * vector;
- }
- #endif
-
- #ifdef VIENNACL_WITH_UBLAS
- // ----------------------------------------------------
- // UBLAS
- //
- template< typename MatrixT, typename VectorT >
- typename viennacl::enable_if< viennacl::is_ublas< typename viennacl::traits::tag_of< MatrixT >::type >::value,
- VectorT>::type
- prod(MatrixT const& matrix, VectorT const& vector)
- {
- // std::cout << "ublas .. " << std::endl;
- return boost::numeric::ublas::prod(matrix, vector);
- }
- #endif
-
-
- // ----------------------------------------------------
- // STL type
- //
-
- // dense matrix-vector product:
- template< typename T, typename A1, typename A2, typename VectorT >
- VectorT
- prod(std::vector< std::vector<T, A1>, A2 > const & matrix, VectorT const& vector)
- {
- VectorT result(matrix.size());
- for (typename std::vector<T, A1>::size_type i=0; i<matrix.size(); ++i)
- {
- result[i] = 0; //we will not assume that VectorT is initialized to zero
- for (typename std::vector<T, A1>::size_type j=0; j<matrix[i].size(); ++j)
- result[i] += matrix[i][j] * vector[j];
- }
- return result;
- }
-
- // sparse matrix-vector product:
- template< typename KEY, typename DATA, typename COMPARE, typename AMAP, typename AVEC, typename VectorT >
- VectorT
- prod(std::vector< std::map<KEY, DATA, COMPARE, AMAP>, AVEC > const& matrix, VectorT const& vector)
- {
- typedef std::vector< std::map<KEY, DATA, COMPARE, AMAP>, AVEC > MatrixType;
-
- VectorT result(matrix.size());
- for (typename MatrixType::size_type i=0; i<matrix.size(); ++i)
- {
- result[i] = 0; //we will not assume that VectorT is initialized to zero
- for (typename std::map<KEY, DATA, COMPARE, AMAP>::const_iterator row_entries = matrix[i].begin();
- row_entries != matrix[i].end();
- ++row_entries)
- result[i] += row_entries->second * vector[row_entries->first];
- }
- return result;
- }
-
-
- /*template< typename MatrixT, typename VectorT >
- VectorT
- prod(MatrixT const& matrix, VectorT const& vector,
- typename viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< MatrixT >::type >::value
- >::type* dummy = 0)
- {
- // std::cout << "std .. " << std::endl;
- return prod_impl(matrix, vector);
- }*/
-
- // ----------------------------------------------------
- // VIENNACL
- //
-
- // standard product:
- template<typename NumericT>
- viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::matrix_base<NumericT>,
- viennacl::op_mat_mat_prod >
- prod(viennacl::matrix_base<NumericT> const & A,
- viennacl::matrix_base<NumericT> const & B)
- {
- return viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::matrix_base<NumericT>,
- viennacl::op_mat_mat_prod >(A, B);
- }
-
- // right factor is a matrix expression:
- template<typename NumericT, typename LhsT, typename RhsT, typename OpT>
- viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- viennacl::op_mat_mat_prod >
- prod(viennacl::matrix_base<NumericT> const & A,
- viennacl::matrix_expression<const LhsT, const RhsT, OpT> const & B)
- {
- return viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- viennacl::op_mat_mat_prod >(A, B);
- }
-
- // left factor is a matrix expression:
- template<typename LhsT, typename RhsT, typename OpT, typename NumericT>
- viennacl::matrix_expression< const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- const viennacl::matrix_base<NumericT>,
- viennacl::op_mat_mat_prod >
- prod(viennacl::matrix_expression<const LhsT, const RhsT, OpT> const & A,
- viennacl::matrix_base<NumericT> const & B)
- {
- return viennacl::matrix_expression< const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- const viennacl::matrix_base<NumericT>,
- viennacl::op_mat_mat_prod >(A, B);
- }
-
-
- // both factors transposed:
- template<typename LhsT1, typename RhsT1, typename OpT1,
- typename LhsT2, typename RhsT2, typename OpT2>
- viennacl::matrix_expression< const viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1>,
- const viennacl::matrix_expression<const LhsT2, const RhsT2, OpT2>,
- viennacl::op_mat_mat_prod >
- prod(viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1> const & A,
- viennacl::matrix_expression<const LhsT2, const RhsT2, OpT2> const & B)
- {
- return viennacl::matrix_expression< const viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1>,
- const viennacl::matrix_expression<const LhsT2, const RhsT2, OpT2>,
- viennacl::op_mat_mat_prod >(A, B);
- }
-
-
-
- // matrix-vector product
- template< typename NumericT>
- viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >
- prod(viennacl::matrix_base<NumericT> const & A,
- viennacl::vector_base<NumericT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >(A, x);
- }
-
- // matrix-vector product (resolve ambiguity)
- template<typename NumericT, typename F>
- viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >
- prod(viennacl::matrix<NumericT, F> const & A,
- viennacl::vector_base<NumericT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >(A, x);
- }
-
- // matrix-vector product (resolve ambiguity)
- template<typename MatrixT, typename NumericT>
- viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >
- prod(viennacl::matrix_range<MatrixT> const & A,
- viennacl::vector_base<NumericT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >(A, x);
- }
-
- // matrix-vector product (resolve ambiguity)
- template<typename MatrixT, typename NumericT>
- viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >
- prod(viennacl::matrix_slice<MatrixT> const & A,
- viennacl::vector_base<NumericT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >(A, x);
- }
-
- // matrix-vector product with matrix expression (including transpose)
- template< typename NumericT, typename LhsT, typename RhsT, typename OpT>
- viennacl::vector_expression< const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >
- prod(viennacl::matrix_expression<const LhsT, const RhsT, OpT> const & A,
- viennacl::vector_base<NumericT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_expression<const LhsT, const RhsT, OpT>,
- const viennacl::vector_base<NumericT>,
- viennacl::op_prod >(A, x);
- }
-
-
- // matrix-vector product with vector expression
- template< typename NumericT, typename LhsT, typename RhsT, typename OpT>
- viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_expression<const LhsT, const RhsT, OpT>,
- viennacl::op_prod >
- prod(viennacl::matrix_base<NumericT> const & A,
- viennacl::vector_expression<const LhsT, const RhsT, OpT> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_base<NumericT>,
- const viennacl::vector_expression<const LhsT, const RhsT, OpT>,
- viennacl::op_prod >(A, x);
- }
-
-
- // matrix-vector product with matrix expression (including transpose) and vector expression
- template<typename LhsT1, typename RhsT1, typename OpT1,
- typename LhsT2, typename RhsT2, typename OpT2>
- viennacl::vector_expression< const viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1>,
- const viennacl::vector_expression<const LhsT2, const RhsT2, OpT2>,
- viennacl::op_prod >
- prod(viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1> const & A,
- viennacl::vector_expression<const LhsT2, const RhsT2, OpT2> const & x)
- {
- return viennacl::vector_expression< const viennacl::matrix_expression<const LhsT1, const RhsT1, OpT1>,
- const viennacl::vector_expression<const LhsT2, const RhsT2, OpT2>,
- viennacl::op_prod >(A, x);
- }
-
-
-
-
- template< typename SparseMatrixType, typename SCALARTYPE>
- typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
- viennacl::matrix_expression<const SparseMatrixType,
- const matrix_base <SCALARTYPE>,
- op_prod >
- >::type
- prod(const SparseMatrixType & sp_mat,
- const viennacl::matrix_base<SCALARTYPE> & d_mat)
- {
- return viennacl::matrix_expression<const SparseMatrixType,
- const viennacl::matrix_base<SCALARTYPE>,
- op_prod >(sp_mat, d_mat);
- }
-
- // right factor is transposed
- template< typename SparseMatrixType, typename SCALARTYPE>
- typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
- viennacl::matrix_expression< const SparseMatrixType,
- const viennacl::matrix_expression<const viennacl::matrix_base<SCALARTYPE>,
- const viennacl::matrix_base<SCALARTYPE>,
- op_trans>,
- viennacl::op_prod >
- >::type
- prod(const SparseMatrixType & A,
- viennacl::matrix_expression<const viennacl::matrix_base<SCALARTYPE>,
- const viennacl::matrix_base<SCALARTYPE>,
- op_trans> const & B)
- {
- return viennacl::matrix_expression< const SparseMatrixType,
- const viennacl::matrix_expression<const viennacl::matrix_base<SCALARTYPE>,
- const viennacl::matrix_base<SCALARTYPE>,
- op_trans>,
- viennacl::op_prod >(A, B);
- }
-
-
- /** @brief Sparse matrix-matrix product with compressed_matrix objects */
- template<typename NumericT>
- viennacl::matrix_expression<const compressed_matrix<NumericT>,
- const compressed_matrix<NumericT>,
- op_prod >
- prod(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> const & B)
- {
- return viennacl::matrix_expression<const compressed_matrix<NumericT>,
- const compressed_matrix<NumericT>,
- op_prod >(A, B);
- }
-
- /** @brief Generic matrix-vector product with user-provided sparse matrix type */
- template<typename SparseMatrixType, typename NumericT>
- vector_expression<const SparseMatrixType,
- const vector_base<NumericT>,
- op_prod >
- prod(const SparseMatrixType & A,
- const vector_base<NumericT> & x)
- {
- return vector_expression<const SparseMatrixType,
- const vector_base<NumericT>,
- op_prod >(A, x);
- }
-
- } // end namespace linalg
-} // end namespace viennacl
-#endif
-
-
-
-
-