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/08 21:40:16 UTC
[19/51] [partial] mahout git commit: (nojira) add native-viennaCL
module to codebase. closes apache/mahout#241
http://git-wip-us.apache.org/repos/asf/mahout/blob/f7c1f802/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/misc_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/misc_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/misc_operations.hpp
new file mode 100644
index 0000000..11061d9
--- /dev/null
+++ b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/misc_operations.hpp
@@ -0,0 +1,80 @@
+#ifndef VIENNACL_LINALG_HOST_BASED_MISC_OPERATIONS_HPP_
+#define VIENNACL_LINALG_HOST_BASED_MISC_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/host_based/misc_operations.hpp
+ @brief Implementations of miscellaneous operations on the CPU using a single thread or OpenMP.
+*/
+
+#include <list>
+
+#include "viennacl/forwards.h"
+#include "viennacl/scalar.hpp"
+#include "viennacl/vector.hpp"
+#include "viennacl/tools/tools.hpp"
+#include "viennacl/linalg/host_based/common.hpp"
+
+namespace viennacl
+{
+namespace linalg
+{
+namespace host_based
+{
+namespace detail
+{
+ template<typename NumericT>
+ void level_scheduling_substitute(vector<NumericT> & vec,
+ viennacl::backend::mem_handle const & row_index_array,
+ viennacl::backend::mem_handle const & row_buffer,
+ viennacl::backend::mem_handle const & col_buffer,
+ viennacl::backend::mem_handle const & element_buffer,
+ vcl_size_t num_rows
+ )
+ {
+ NumericT * vec_buf = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(vec.handle());
+
+ unsigned int const * elim_row_index = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(row_index_array);
+ unsigned int const * elim_row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(row_buffer);
+ unsigned int const * elim_col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(col_buffer);
+ NumericT const * elim_elements = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(element_buffer);
+
+#ifdef VIENNACL_WITH_OPENMP
+ #pragma omp parallel for
+#endif
+ for (long row=0; row < static_cast<long>(num_rows); ++row)
+ {
+ unsigned int eq_row = elim_row_index[row];
+ unsigned int row_end = elim_row_buffer[row+1];
+ NumericT vec_entry = vec_buf[eq_row];
+
+ for (vcl_size_t j = elim_row_buffer[row]; j < row_end; ++j)
+ vec_entry -= vec_buf[elim_col_buffer[j]] * elim_elements[j];
+
+ vec_buf[eq_row] = vec_entry;
+ }
+
+ }
+}
+
+} // namespace host_based
+} //namespace linalg
+} //namespace viennacl
+
+
+#endif
http://git-wip-us.apache.org/repos/asf/mahout/blob/f7c1f802/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/nmf_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/nmf_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/nmf_operations.hpp
new file mode 100644
index 0000000..bb6557f
--- /dev/null
+++ b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/nmf_operations.hpp
@@ -0,0 +1,247 @@
+#ifndef VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_
+#define VIENNACL_LINALG_HOST_BASED_NMF_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/host_based/vector_operations.hpp
+ @brief Implementations of NMF operations using a plain single-threaded or OpenMP-enabled execution on CPU
+ */
+
+#include "viennacl/vector.hpp"
+#include "viennacl/matrix.hpp"
+#include "viennacl/linalg/prod.hpp"
+#include "viennacl/linalg/norm_2.hpp"
+#include "viennacl/linalg/norm_frobenius.hpp"
+
+#include "viennacl/linalg/host_based/common.hpp"
+
+namespace viennacl
+{
+namespace linalg
+{
+
+/** @brief Configuration class for the nonnegative-matrix-factorization algorithm. Specify tolerances, maximum iteration counts, etc., here. */
+class nmf_config
+{
+public:
+ nmf_config(double val_epsilon = 1e-4, double val_epsilon_stagnation = 1e-5,
+ vcl_size_t num_max_iters = 10000, vcl_size_t num_check_iters = 100) :
+ eps_(val_epsilon), stagnation_eps_(val_epsilon_stagnation), max_iters_(num_max_iters), check_after_steps_(
+ (num_check_iters > 0) ? num_check_iters : 1), print_relative_error_(false), iters_(0)
+ {
+ }
+
+ /** @brief Returns the relative tolerance for convergence */
+ double tolerance() const
+ {
+ return eps_;
+ }
+
+ /** @brief Sets the relative tolerance for convergence, i.e. norm(V - W * H) / norm(V - W_init * H_init) */
+ void tolerance(double e)
+ {
+ eps_ = e;
+ }
+
+ /** @brief Relative tolerance for the stagnation check */
+ double stagnation_tolerance() const
+ {
+ return stagnation_eps_;
+ }
+
+ /** @brief Sets the tolerance for the stagnation check (i.e. the minimum required relative change of the residual between two iterations) */
+ void stagnation_tolerance(double e)
+ {
+ stagnation_eps_ = e;
+ }
+
+ /** @brief Returns the maximum number of iterations for the NMF algorithm */
+ vcl_size_t max_iterations() const
+ {
+ return max_iters_;
+ }
+ /** @brief Sets the maximum number of iterations for the NMF algorithm */
+ void max_iterations(vcl_size_t m)
+ {
+ max_iters_ = m;
+ }
+
+ /** @brief Returns the number of iterations of the last NMF run using this configuration object */
+ vcl_size_t iters() const
+ {
+ return iters_;
+ }
+
+ /** @brief Number of steps after which the convergence of NMF should be checked (again) */
+ vcl_size_t check_after_steps() const
+ {
+ return check_after_steps_;
+ }
+ /** @brief Set the number of steps after which the convergence of NMF should be checked (again) */
+ void check_after_steps(vcl_size_t c)
+ {
+ if (c > 0)
+ check_after_steps_ = c;
+ }
+
+ /** @brief Returns the flag specifying whether the relative tolerance should be printed in each iteration */
+ bool print_relative_error() const
+ {
+ return print_relative_error_;
+ }
+ /** @brief Specify whether the relative error should be printed at each convergence check after 'num_check_iters' steps */
+ void print_relative_error(bool b)
+ {
+ print_relative_error_ = b;
+ }
+
+ template<typename ScalarType>
+ friend void nmf(viennacl::matrix_base<ScalarType> const & V,
+ viennacl::matrix_base<ScalarType> & W, viennacl::matrix_base<ScalarType> & H,
+ nmf_config const & conf);
+
+private:
+ double eps_;
+ double stagnation_eps_;
+ vcl_size_t max_iters_;
+ vcl_size_t check_after_steps_;
+ bool print_relative_error_;
+public:
+ mutable vcl_size_t iters_;
+};
+
+namespace host_based
+{
+ /** @brief Missing OpenMP kernel for nonnegative matrix factorization of a dense matrices. */
+ template<typename NumericT>
+ void el_wise_mul_div(NumericT * matrix1,
+ NumericT const * matrix2,
+ NumericT const * matrix3, vcl_size_t size)
+ {
+#ifdef VIENNACL_WITH_OPENMP
+#pragma omp parallel for
+#endif
+ for (long i2 = 0; i2 < long(size); i2++)
+ {
+ vcl_size_t i = vcl_size_t(i2);
+ NumericT val = matrix1[i] * matrix2[i];
+ NumericT divisor = matrix3[i];
+ matrix1[i] = (divisor > (NumericT) 0.00001) ? (val / divisor) : (NumericT) 0;
+ }
+ }
+
+ /** @brief The nonnegative matrix factorization (approximation) algorithm as suggested by Lee and Seung. Factorizes a matrix V with nonnegative entries into matrices W and H such that ||V - W*H|| is minimized.
+ *
+ * @param V Input matrix
+ * @param W First factor
+ * @param H Second factor
+ * @param conf A configuration object holding tolerances and the like
+ */
+ template<typename NumericT>
+ void nmf(viennacl::matrix_base<NumericT> const & V,
+ viennacl::matrix_base<NumericT> & W,
+ viennacl::matrix_base<NumericT> & H,
+ viennacl::linalg::nmf_config const & conf)
+ {
+ vcl_size_t k = W.size2();
+ conf.iters_ = 0;
+
+ if (viennacl::linalg::norm_frobenius(W) <= 0)
+ W = viennacl::scalar_matrix<NumericT>(W.size1(), W.size2(), NumericT(1.0));
+
+ if (viennacl::linalg::norm_frobenius(H) <= 0)
+ H = viennacl::scalar_matrix<NumericT>(H.size1(), H.size2(), NumericT(1.0));
+
+ viennacl::matrix_base<NumericT> wn(V.size1(), k, W.row_major());
+ viennacl::matrix_base<NumericT> wd(V.size1(), k, W.row_major());
+ viennacl::matrix_base<NumericT> wtmp(V.size1(), V.size2(), W.row_major());
+
+ viennacl::matrix_base<NumericT> hn(k, V.size2(), H.row_major());
+ viennacl::matrix_base<NumericT> hd(k, V.size2(), H.row_major());
+ viennacl::matrix_base<NumericT> htmp(k, k, H.row_major());
+
+ viennacl::matrix_base<NumericT> appr(V.size1(), V.size2(), V.row_major());
+
+ NumericT last_diff = 0;
+ NumericT diff_init = 0;
+ bool stagnation_flag = false;
+
+ for (vcl_size_t i = 0; i < conf.max_iterations(); i++)
+ {
+ conf.iters_ = i + 1;
+
+ hn = viennacl::linalg::prod(trans(W), V);
+ htmp = viennacl::linalg::prod(trans(W), W);
+ hd = viennacl::linalg::prod(htmp, H);
+
+ NumericT * data_H = detail::extract_raw_pointer<NumericT>(H);
+ NumericT * data_hn = detail::extract_raw_pointer<NumericT>(hn);
+ NumericT * data_hd = detail::extract_raw_pointer<NumericT>(hd);
+
+ viennacl::linalg::host_based::el_wise_mul_div(data_H, data_hn, data_hd, H.internal_size1() * H.internal_size2());
+
+ wn = viennacl::linalg::prod(V, trans(H));
+ wtmp = viennacl::linalg::prod(W, H);
+ wd = viennacl::linalg::prod(wtmp, trans(H));
+
+ NumericT * data_W = detail::extract_raw_pointer<NumericT>(W);
+ NumericT * data_wn = detail::extract_raw_pointer<NumericT>(wn);
+ NumericT * data_wd = detail::extract_raw_pointer<NumericT>(wd);
+
+ viennacl::linalg::host_based::el_wise_mul_div(data_W, data_wn, data_wd, W.internal_size1() * W.internal_size2());
+
+ if (i % conf.check_after_steps() == 0) //check for convergence
+ {
+ appr = viennacl::linalg::prod(W, H);
+
+ appr -= V;
+ NumericT diff_val = viennacl::linalg::norm_frobenius(appr);
+
+ if (i == 0)
+ diff_init = diff_val;
+
+ if (conf.print_relative_error())
+ std::cout << diff_val / diff_init << std::endl;
+
+ // Approximation check
+ if (diff_val / diff_init < conf.tolerance())
+ break;
+
+ // Stagnation check
+ if (std::fabs(diff_val - last_diff) / (diff_val * NumericT(conf.check_after_steps())) < conf.stagnation_tolerance()) //avoid situations where convergence stagnates
+ {
+ if (stagnation_flag) // iteration stagnates (two iterates with no notable progress)
+ break;
+ else
+ // record stagnation in this iteration
+ stagnation_flag = true;
+ } else
+ // good progress in this iteration, so unset stagnation flag
+ stagnation_flag = false;
+
+ // prepare for next iterate:
+ last_diff = diff_val;
+ }
+ }
+ }
+
+} //namespace host_based
+} //namespace linalg
+} //namespace viennacl
+
+#endif /* VIENNACL_LINALG_HOST_BASED_NMF_OPERATIONS_HPP_ */
http://git-wip-us.apache.org/repos/asf/mahout/blob/f7c1f802/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/scalar_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/scalar_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/scalar_operations.hpp
new file mode 100644
index 0000000..f8a1f3b
--- /dev/null
+++ b/native-viennaCL/src/main/cpp/viennacl/linalg/host_based/scalar_operations.hpp
@@ -0,0 +1,162 @@
+#ifndef VIENNACL_LINALG_HOST_BASED_SCALAR_OPERATIONS_HPP_
+#define VIENNACL_LINALG_HOST_BASED_SCALAR_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/host_based/scalar_operations.hpp
+ @brief Implementations of scalar operations using a plain single-threaded or OpenMP-enabled execution on CPU
+*/
+
+#include "viennacl/forwards.h"
+#include "viennacl/tools/tools.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/stride.hpp"
+#include "viennacl/linalg/host_based/common.hpp"
+
+namespace viennacl
+{
+namespace linalg
+{
+namespace host_based
+{
+template<typename ScalarT1,
+ typename ScalarT2, typename FactorT>
+typename viennacl::enable_if< viennacl::is_scalar<ScalarT1>::value
+ && viennacl::is_scalar<ScalarT2>::value
+ && viennacl::is_any_scalar<FactorT>::value
+ >::type
+as(ScalarT1 & s1,
+ ScalarT2 const & s2, FactorT const & alpha, vcl_size_t /*len_alpha*/, bool reciprocal_alpha, bool flip_sign_alpha)
+{
+ typedef typename viennacl::result_of::cpu_value_type<ScalarT1>::type value_type;
+
+ value_type * data_s1 = detail::extract_raw_pointer<value_type>(s1);
+ value_type const * data_s2 = detail::extract_raw_pointer<value_type>(s2);
+
+ value_type data_alpha = alpha;
+ if (flip_sign_alpha)
+ data_alpha = -data_alpha;
+ if (reciprocal_alpha)
+ data_alpha = static_cast<value_type>(1) / data_alpha;
+
+ *data_s1 = *data_s2 * data_alpha;
+}
+
+
+template<typename ScalarT1,
+ typename ScalarT2, typename FactorT2,
+ typename ScalarT3, typename FactorT3>
+typename viennacl::enable_if< viennacl::is_scalar<ScalarT1>::value
+ && viennacl::is_scalar<ScalarT2>::value
+ && viennacl::is_scalar<ScalarT3>::value
+ && viennacl::is_any_scalar<FactorT2>::value
+ && viennacl::is_any_scalar<FactorT3>::value
+ >::type
+asbs(ScalarT1 & s1,
+ ScalarT2 const & s2, FactorT2 const & alpha, vcl_size_t /*len_alpha*/, bool reciprocal_alpha, bool flip_sign_alpha,
+ ScalarT3 const & s3, FactorT3 const & beta, vcl_size_t /*len_beta*/, bool reciprocal_beta, bool flip_sign_beta)
+{
+ typedef typename viennacl::result_of::cpu_value_type<ScalarT1>::type value_type;
+
+ value_type * data_s1 = detail::extract_raw_pointer<value_type>(s1);
+ value_type const * data_s2 = detail::extract_raw_pointer<value_type>(s2);
+ value_type const * data_s3 = detail::extract_raw_pointer<value_type>(s3);
+
+ value_type data_alpha = alpha;
+ if (flip_sign_alpha)
+ data_alpha = -data_alpha;
+ if (reciprocal_alpha)
+ data_alpha = static_cast<value_type>(1) / data_alpha;
+
+ value_type data_beta = beta;
+ if (flip_sign_beta)
+ data_beta = -data_beta;
+ if (reciprocal_beta)
+ data_beta = static_cast<value_type>(1) / data_beta;
+
+ *data_s1 = *data_s2 * data_alpha + *data_s3 * data_beta;
+}
+
+
+template<typename ScalarT1,
+ typename ScalarT2, typename FactorT2,
+ typename ScalarT3, typename FactorT3>
+typename viennacl::enable_if< viennacl::is_scalar<ScalarT1>::value
+ && viennacl::is_scalar<ScalarT2>::value
+ && viennacl::is_scalar<ScalarT3>::value
+ && viennacl::is_any_scalar<FactorT2>::value
+ && viennacl::is_any_scalar<FactorT3>::value
+ >::type
+asbs_s(ScalarT1 & s1,
+ ScalarT2 const & s2, FactorT2 const & alpha, vcl_size_t /*len_alpha*/, bool reciprocal_alpha, bool flip_sign_alpha,
+ ScalarT3 const & s3, FactorT3 const & beta, vcl_size_t /*len_beta*/, bool reciprocal_beta, bool flip_sign_beta)
+{
+ typedef typename viennacl::result_of::cpu_value_type<ScalarT1>::type value_type;
+
+ value_type * data_s1 = detail::extract_raw_pointer<value_type>(s1);
+ value_type const * data_s2 = detail::extract_raw_pointer<value_type>(s2);
+ value_type const * data_s3 = detail::extract_raw_pointer<value_type>(s3);
+
+ value_type data_alpha = alpha;
+ if (flip_sign_alpha)
+ data_alpha = -data_alpha;
+ if (reciprocal_alpha)
+ data_alpha = static_cast<value_type>(1) / data_alpha;
+
+ value_type data_beta = beta;
+ if (flip_sign_beta)
+ data_beta = -data_beta;
+ if (reciprocal_beta)
+ data_beta = static_cast<value_type>(1) / data_beta;
+
+ *data_s1 += *data_s2 * data_alpha + *data_s3 * data_beta;
+}
+
+
+/** @brief Swaps the contents of two scalars, data is copied
+*
+* @param s1 The first scalar
+* @param s2 The second scalar
+*/
+template<typename ScalarT1, typename ScalarT2>
+typename viennacl::enable_if< viennacl::is_scalar<ScalarT1>::value
+ && viennacl::is_scalar<ScalarT2>::value
+ >::type
+swap(ScalarT1 & s1, ScalarT2 & s2)
+{
+ typedef typename viennacl::result_of::cpu_value_type<ScalarT1>::type value_type;
+
+ value_type * data_s1 = detail::extract_raw_pointer<value_type>(s1);
+ value_type * data_s2 = detail::extract_raw_pointer<value_type>(s2);
+
+ value_type temp = *data_s2;
+ *data_s2 = *data_s1;
+ *data_s1 = temp;
+}
+
+
+
+} //namespace host_based
+} //namespace linalg
+} //namespace viennacl
+
+
+#endif