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:21 UTC
[16/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/ilu.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/ilu.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/ilu.hpp
deleted file mode 100644
index 540ff82..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/ilu.hpp
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef VIENNACL_LINALG_ILU_HPP_
-#define VIENNACL_LINALG_ILU_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/ilu.hpp
- @brief Implementations of incomplete factorization preconditioners. Convenience header file.
-*/
-
-#include "viennacl/linalg/detail/ilu/ilut.hpp"
-#include "viennacl/linalg/detail/ilu/ilu0.hpp"
-#include "viennacl/linalg/detail/ilu/block_ilu.hpp"
-#include "viennacl/linalg/detail/ilu/chow_patel_ilu.hpp"
-
-#endif
-
-
-
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/ilu_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/ilu_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/ilu_operations.hpp
deleted file mode 100644
index febd347..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/ilu_operations.hpp
+++ /dev/null
@@ -1,334 +0,0 @@
-#ifndef VIENNACL_LINALG_ILU_OPERATIONS_HPP_
-#define VIENNACL_LINALG_ILU_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 PDF manual)
-
- License: MIT (X11), see file LICENSE in the base directory
-============================================================================= */
-
-/** @file viennacl/linalg/ilu_operations.hpp
- @brief Implementations of specialized routines for the Chow-Patel parallel ILU preconditioner
-*/
-
-#include "viennacl/forwards.h"
-#include "viennacl/range.hpp"
-#include "viennacl/scalar.hpp"
-#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/handle.hpp"
-#include "viennacl/traits/stride.hpp"
-#include "viennacl/linalg/host_based/ilu_operations.hpp"
-
-#ifdef VIENNACL_WITH_OPENCL
- #include "viennacl/linalg/opencl/ilu_operations.hpp"
-#endif
-
-#ifdef VIENNACL_WITH_CUDA
- #include "viennacl/linalg/cuda/ilu_operations.hpp"
-#endif
-
-namespace viennacl
-{
-namespace linalg
-{
-
-/** @brief Extracts the lower triangular part L from A.
- *
- * Diagonal of L is stored explicitly in order to enable better code reuse.
- *
- */
-template<typename NumericT>
-void extract_L(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> & L)
-{
- switch (viennacl::traits::handle(A).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::extract_L(A, L);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::extract_L(A, L);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::extract_L(A, L);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L accordingly.
- *
- * Since A should not be modified (const-correctness), updates are in L.
- *
- */
-template<typename NumericT>
-void icc_scale(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> & L)
-{
- switch (viennacl::traits::handle(A).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::icc_scale(A, L);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::icc_scale(A, L);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::icc_scale(A, L);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Performs one nonlinear relaxation step in the Chow-Patel-ICC (cf. Algorithm 3 in paper, but for L rather than U)
- *
- * We use a fully synchronous (Jacobi-like) variant, because asynchronous methods as described in the paper are a nightmare to debug
- * (and particularly funny if they sometimes fail, sometimes not)
- *
- * @param L Factor L to be updated for the incomplete Cholesky factorization
- * @param aij_L Lower triangular potion from system matrix
- */
-template<typename NumericT>
-void icc_chow_patel_sweep(compressed_matrix<NumericT> & L,
- vector<NumericT> & aij_L)
-{
- switch (viennacl::traits::handle(L).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::icc_chow_patel_sweep(L, aij_L);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::icc_chow_patel_sweep(L, aij_L);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::icc_chow_patel_sweep(L, aij_L);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-
-//////////////////////// ILU ////////////////////
-
-/** @brief Extracts the lower triangular part L and the upper triangular part U from A.
- *
- * Diagonals of L and U are stored explicitly in order to enable better code reuse.
- *
- */
-template<typename NumericT>
-void extract_LU(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> & L,
- compressed_matrix<NumericT> & U)
-{
- switch (viennacl::traits::handle(A).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::extract_LU(A, L, U);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::extract_LU(A, L, U);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::extract_LU(A, L, U);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L and U accordingly.
- *
- * Since A should not be modified (const-correctness), updates are in L and U.
- *
- */
-template<typename NumericT>
-void ilu_scale(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> & L,
- compressed_matrix<NumericT> & U)
-{
- switch (viennacl::traits::handle(A).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::ilu_scale(A, L, U);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::ilu_scale(A, L, U);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::ilu_scale(A, L, U);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Transposition B <- A^T, where the aij-vector is permuted in the same way as the value array in A when assigned to B
- *
- * @param A Input matrix to be transposed
- * @param B Output matrix containing the transposed matrix
- */
-template<typename NumericT>
-void ilu_transpose(compressed_matrix<NumericT> const & A,
- compressed_matrix<NumericT> & B)
-{
- viennacl::context orig_ctx = viennacl::traits::context(A);
- viennacl::context cpu_ctx(viennacl::MAIN_MEMORY);
- (void)orig_ctx;
- (void)cpu_ctx;
-
- viennacl::compressed_matrix<NumericT> A_host(0, 0, 0, cpu_ctx);
- (void)A_host;
-
- switch (viennacl::traits::handle(A).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::ilu_transpose(A, B);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- A_host = A;
- B.switch_memory_context(cpu_ctx);
- viennacl::linalg::host_based::ilu_transpose(A_host, B);
- B.switch_memory_context(orig_ctx);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- A_host = A;
- B.switch_memory_context(cpu_ctx);
- viennacl::linalg::host_based::ilu_transpose(A_host, B);
- B.switch_memory_context(orig_ctx);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-
-/** @brief Performs one nonlinear relaxation step in the Chow-Patel-ILU (cf. Algorithm 2 in paper)
- *
- * We use a fully synchronous (Jacobi-like) variant, because asynchronous methods as described in the paper are a nightmare to debug
- * (and particularly funny if they sometimes fail, sometimes not)
- *
- * @param L Lower-triangular matrix L in LU factorization
- * @param aij_L Lower-triangular matrix L from A
- * @param U_trans Upper-triangular matrix U in CSC-storage, which is the same as U^trans in CSR-storage
- * @param aij_U_trans Upper-triangular matrix from A in CSC-storage, which is the same as U^trans in CSR-storage
- */
-template<typename NumericT>
-void ilu_chow_patel_sweep(compressed_matrix<NumericT> & L,
- vector<NumericT> const & aij_L,
- compressed_matrix<NumericT> & U_trans,
- vector<NumericT> const & aij_U_trans)
-{
- switch (viennacl::traits::handle(L).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::ilu_chow_patel_sweep(L, aij_L, U_trans, aij_U_trans);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::ilu_chow_patel_sweep(L, aij_L, U_trans, aij_U_trans);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::ilu_chow_patel_sweep(L, aij_L, U_trans, aij_U_trans);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Extracts the lower triangular part L and the upper triangular part U from A.
- *
- * Diagonals of L and U are stored explicitly in order to enable better code reuse.
- *
- */
-template<typename NumericT>
-void ilu_form_neumann_matrix(compressed_matrix<NumericT> & R,
- vector<NumericT> & diag_R)
-{
- switch (viennacl::traits::handle(R).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::ilu_form_neumann_matrix(R, diag_R);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::ilu_form_neumann_matrix(R, diag_R);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::ilu_form_neumann_matrix(R, diag_R);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-} //namespace linalg
-} //namespace viennacl
-
-
-#endif
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/inner_prod.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/inner_prod.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/inner_prod.hpp
deleted file mode 100644
index b31a82a..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/inner_prod.hpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#ifndef VIENNACL_LINALG_INNER_PROD_HPP_
-#define VIENNACL_LINALG_INNER_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/inner_prod.hpp
- @brief Generic interface for the computation of inner products. See viennacl/linalg/vector_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 "viennacl/meta/result_of.hpp"
-
-namespace viennacl
-{
-//
-// generic inner_prod function
-// uses tag dispatch to identify which algorithm
-// should be called
-//
-namespace linalg
-{
-
-#ifdef VIENNACL_WITH_ARMADILLO
-// ----------------------------------------------------
-// Armadillo
-//
-template<typename NumericT>
-NumericT inner_prod(arma::Col<NumericT> const& v1, arma::Col<NumericT> const& v2)
-{
- return dot(v1, v2);
-}
-#endif
-
-#ifdef VIENNACL_WITH_EIGEN
-// ----------------------------------------------------
-// EIGEN
-//
-template<typename VectorT1, typename VectorT2>
-typename viennacl::enable_if< viennacl::is_eigen< typename viennacl::traits::tag_of< VectorT1 >::type >::value,
- typename VectorT1::RealScalar>::type
-inner_prod(VectorT1 const & v1, VectorT2 const & v2)
-{
- //std::cout << "eigen .. " << std::endl;
- return v1.dot(v2);
-}
-#endif
-
-#ifdef VIENNACL_WITH_MTL4
-// ----------------------------------------------------
-// MTL4
-//
-template<typename VectorT1, typename VectorT2>
-typename viennacl::enable_if< viennacl::is_mtl4< typename viennacl::traits::tag_of< VectorT1 >::type >::value,
- typename VectorT1::value_type>::type
-inner_prod(VectorT1 const & v1, VectorT2 const & v2)
-{
- //std::cout << "mtl4 .. " << std::endl;
- return mtl::dot(v1, v2);
-}
-#endif
-
-#ifdef VIENNACL_WITH_UBLAS
-// ----------------------------------------------------
-// UBLAS
-//
-template<typename VectorT1, typename VectorT2>
-typename viennacl::enable_if< viennacl::is_ublas< typename viennacl::traits::tag_of< VectorT1 >::type >::value,
- typename VectorT1::value_type>::type
-inner_prod(VectorT1 const & v1, VectorT2 const & v2)
-{
- //std::cout << "ublas .. " << std::endl;
- return boost::numeric::ublas::inner_prod(v1, v2);
-}
-#endif
-
-// ----------------------------------------------------
-// STL
-//
-template<typename VectorT1, typename VectorT2>
-typename viennacl::enable_if< viennacl::is_stl< typename viennacl::traits::tag_of< VectorT1 >::type >::value,
- typename VectorT1::value_type>::type
-inner_prod(VectorT1 const & v1, VectorT2 const & v2)
-{
- assert(v1.size() == v2.size() && bool("Vector sizes mismatch"));
- //std::cout << "stl .. " << std::endl;
- typename VectorT1::value_type result = 0;
- for (typename VectorT1::size_type i=0; i<v1.size(); ++i)
- result += v1[i] * v2[i];
-
- return result;
-}
-
-// ----------------------------------------------------
-// VIENNACL
-//
-template<typename NumericT>
-viennacl::scalar_expression< const vector_base<NumericT>, const vector_base<NumericT>, viennacl::op_inner_prod >
-inner_prod(vector_base<NumericT> const & vector1,
- vector_base<NumericT> const & vector2)
-{
- //std::cout << "viennacl .. " << std::endl;
- return viennacl::scalar_expression< const vector_base<NumericT>,
- const vector_base<NumericT>,
- viennacl::op_inner_prod >(vector1, vector2);
-}
-
-
-// expression on lhs:
-template< typename LHS, typename RHS, typename OP, typename NumericT>
-viennacl::scalar_expression< const viennacl::vector_expression<LHS, RHS, OP>,
- const vector_base<NumericT>,
- viennacl::op_inner_prod >
-inner_prod(viennacl::vector_expression<LHS, RHS, OP> const & vector1,
- vector_base<NumericT> const & vector2)
-{
- //std::cout << "viennacl .. " << std::endl;
- return viennacl::scalar_expression< const viennacl::vector_expression<LHS, RHS, OP>,
- const vector_base<NumericT>,
- viennacl::op_inner_prod >(vector1, vector2);
-}
-
-// expression on rhs:
-template<typename NumericT, typename LHS, typename RHS, typename OP>
-viennacl::scalar_expression< const vector_base<NumericT>,
- const viennacl::vector_expression<LHS, RHS, OP>,
- viennacl::op_inner_prod >
-inner_prod(vector_base<NumericT> const & vector1,
- viennacl::vector_expression<LHS, RHS, OP> const & vector2)
-{
- //std::cout << "viennacl .. " << std::endl;
- return viennacl::scalar_expression< const vector_base<NumericT>,
- const viennacl::vector_expression<LHS, RHS, OP>,
- viennacl::op_inner_prod >(vector1, vector2);
-}
-
-// expression on lhs and rhs:
-template<typename LHS1, typename RHS1, typename OP1,
- typename LHS2, typename RHS2, typename OP2>
-viennacl::scalar_expression< const viennacl::vector_expression<LHS1, RHS1, OP1>,
- const viennacl::vector_expression<LHS2, RHS2, OP2>,
- viennacl::op_inner_prod >
-inner_prod(viennacl::vector_expression<LHS1, RHS1, OP1> const & vector1,
- viennacl::vector_expression<LHS2, RHS2, OP2> const & vector2)
-{
- //std::cout << "viennacl .. " << std::endl;
- return viennacl::scalar_expression< const viennacl::vector_expression<LHS1, RHS1, OP1>,
- const viennacl::vector_expression<LHS2, RHS2, OP2>,
- viennacl::op_inner_prod >(vector1, vector2);
-}
-
-
-// Multiple inner products:
-template<typename NumericT>
-viennacl::vector_expression< const vector_base<NumericT>, const vector_tuple<NumericT>, viennacl::op_inner_prod >
-inner_prod(vector_base<NumericT> const & x,
- vector_tuple<NumericT> const & y_tuple)
-{
- return viennacl::vector_expression< const vector_base<NumericT>,
- const vector_tuple<NumericT>,
- viennacl::op_inner_prod >(x, y_tuple);
-}
-
-
-} // 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/iterative_operations.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/iterative_operations.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/iterative_operations.hpp
deleted file mode 100644
index 78a813d..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/iterative_operations.hpp
+++ /dev/null
@@ -1,425 +0,0 @@
-#ifndef VIENNACL_LINALG_ITERATIVE_OPERATIONS_HPP_
-#define VIENNACL_LINALG_ITERATIVE_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/iterative_operations.hpp
- @brief Implementations of specialized routines for the iterative solvers.
-*/
-
-#include "viennacl/forwards.h"
-#include "viennacl/range.hpp"
-#include "viennacl/scalar.hpp"
-#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/handle.hpp"
-#include "viennacl/traits/stride.hpp"
-#include "viennacl/linalg/host_based/iterative_operations.hpp"
-
-#ifdef VIENNACL_WITH_OPENCL
- #include "viennacl/linalg/opencl/iterative_operations.hpp"
-#endif
-
-#ifdef VIENNACL_WITH_CUDA
- #include "viennacl/linalg/cuda/iterative_operations.hpp"
-#endif
-
-namespace viennacl
-{
-namespace linalg
-{
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined CG algorithm.
- *
- * This routines computes for vectors 'result', 'p', 'r', 'Ap':
- * result += alpha * p;
- * r -= alpha * Ap;
- * p = r + beta * p;
- * and runs the parallel reduction stage for computing inner_prod(r,r)
- */
-template<typename NumericT>
-void pipelined_cg_vector_update(vector_base<NumericT> & result,
- NumericT alpha,
- vector_base<NumericT> & p,
- vector_base<NumericT> & r,
- vector_base<NumericT> const & Ap,
- NumericT beta,
- vector_base<NumericT> & inner_prod_buffer)
-{
- switch (viennacl::traits::handle(result).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_cg_vector_update(result, alpha, p, r, Ap, beta, inner_prod_buffer);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_cg_vector_update(result, alpha, p, r, Ap, beta, inner_prod_buffer);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_cg_vector_update(result, alpha, p, r, Ap, beta, inner_prod_buffer);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined CG algorithm.
- *
- * This routines computes for a matrix A and vectors 'p' and 'Ap':
- * Ap = prod(A, p);
- * and computes the two reduction stages for computing inner_prod(p,Ap), inner_prod(Ap,Ap)
- */
-template<typename MatrixT, typename NumericT>
-void pipelined_cg_prod(MatrixT const & A,
- vector_base<NumericT> const & p,
- vector_base<NumericT> & Ap,
- vector_base<NumericT> & inner_prod_buffer)
-{
- switch (viennacl::traits::handle(p).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_cg_prod(A, p, Ap, inner_prod_buffer);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_cg_prod(A, p, Ap, inner_prod_buffer);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_cg_prod(A, p, Ap, inner_prod_buffer);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-////////////////////////////////////////////
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined CG algorithm.
- *
- * This routines computes for vectors 's', 'r', 'Ap':
- * s = r - alpha * Ap
- * with alpha obtained from a reduction step on the 0th and the 3rd out of 6 chunks in inner_prod_buffer
- * and runs the parallel reduction stage for computing inner_prod(s,s)
- */
-template<typename NumericT>
-void pipelined_bicgstab_update_s(vector_base<NumericT> & s,
- vector_base<NumericT> & r,
- vector_base<NumericT> const & Ap,
- vector_base<NumericT> & inner_prod_buffer,
- vcl_size_t buffer_chunk_size,
- vcl_size_t buffer_chunk_offset)
-{
- switch (viennacl::traits::handle(s).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_bicgstab_update_s(s, r, Ap, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_bicgstab_update_s(s, r, Ap, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_bicgstab_update_s(s, r, Ap, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined BiCGStab algorithm.
- *
- * x_{j+1} = x_j + alpha * p_j + omega * s_j
- * r_{j+1} = s_j - omega * t_j
- * p_{j+1} = r_{j+1} + beta * (p_j - omega * q_j)
- * and compute first stage of r_dot_r0 = <r_{j+1}, r_o^*> for use in next iteration
- */
-template<typename NumericT>
-void pipelined_bicgstab_vector_update(vector_base<NumericT> & result, NumericT alpha, vector_base<NumericT> & p, NumericT omega, vector_base<NumericT> const & s,
- vector_base<NumericT> & residual, vector_base<NumericT> const & As,
- NumericT beta, vector_base<NumericT> const & Ap,
- vector_base<NumericT> const & r0star,
- vector_base<NumericT> & inner_prod_buffer,
- vcl_size_t buffer_chunk_size)
-{
- switch (viennacl::traits::handle(s).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_bicgstab_vector_update(result, alpha, p, omega, s, residual, As, beta, Ap, r0star, inner_prod_buffer, buffer_chunk_size);
- break;
- #ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_bicgstab_vector_update(result, alpha, p, omega, s, residual, As, beta, Ap, r0star, inner_prod_buffer, buffer_chunk_size);
- break;
- #endif
- #ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_bicgstab_vector_update(result, alpha, p, omega, s, residual, As, beta, Ap, r0star, inner_prod_buffer, buffer_chunk_size);
- break;
- #endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined CG algorithm.
- *
- * This routines computes for a matrix A and vectors 'p' and 'Ap':
- * Ap = prod(A, p);
- * and computes the two reduction stages for computing inner_prod(p,Ap), inner_prod(Ap,Ap)
- */
-template<typename MatrixT, typename NumericT>
-void pipelined_bicgstab_prod(MatrixT const & A,
- vector_base<NumericT> const & p,
- vector_base<NumericT> & Ap,
- vector_base<NumericT> const & r0star,
- vector_base<NumericT> & inner_prod_buffer,
- vcl_size_t buffer_chunk_size,
- vcl_size_t buffer_chunk_offset)
-{
- switch (viennacl::traits::handle(p).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_bicgstab_prod(A, p, Ap, r0star, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_bicgstab_prod(A, p, Ap, r0star, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_bicgstab_prod(A, p, Ap, r0star, inner_prod_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-////////////////////////////////////////////
-
-/** @brief Performs a vector normalization needed for an efficient pipelined GMRES algorithm.
- *
- * This routines computes for vectors 'r', 'v_k':
- * Second reduction step for ||v_k||
- * v_k /= ||v_k||
- * First reduction step for <r, v_k>
- */
-template <typename T>
-void pipelined_gmres_normalize_vk(vector_base<T> & v_k,
- vector_base<T> const & residual,
- vector_base<T> & R_buffer,
- vcl_size_t offset_in_R,
- vector_base<T> const & inner_prod_buffer,
- vector_base<T> & r_dot_vk_buffer,
- vcl_size_t buffer_chunk_size,
- vcl_size_t buffer_chunk_offset)
-{
- switch (viennacl::traits::handle(v_k).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_gmres_normalize_vk(v_k, residual, R_buffer, offset_in_R, inner_prod_buffer, r_dot_vk_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_gmres_normalize_vk(v_k, residual, R_buffer, offset_in_R, inner_prod_buffer, r_dot_vk_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_gmres_normalize_vk(v_k, residual, R_buffer, offset_in_R, inner_prod_buffer, r_dot_vk_buffer, buffer_chunk_size, buffer_chunk_offset);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-
-/** @brief Computes the first reduction stage for multiple inner products <v_i, v_k>, i=0..k-1
- *
- * All vectors v_i are stored column-major in the array 'device_krylov_basis', where each vector has an actual length 'v_k_size', but might be padded to have 'v_k_internal_size'
- */
-template <typename T>
-void pipelined_gmres_gram_schmidt_stage1(vector_base<T> const & device_krylov_basis,
- vcl_size_t v_k_size,
- vcl_size_t v_k_internal_size,
- vcl_size_t k,
- vector_base<T> & vi_in_vk_buffer,
- vcl_size_t buffer_chunk_size)
-{
- switch (viennacl::traits::handle(device_krylov_basis).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_gmres_gram_schmidt_stage1(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, buffer_chunk_size);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_gmres_gram_schmidt_stage1(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, buffer_chunk_size);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_gmres_gram_schmidt_stage1(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, buffer_chunk_size);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-/** @brief Computes the second reduction stage for multiple inner products <v_i, v_k>, i=0..k-1, then updates v_k -= <v_i, v_k> v_i and computes the first reduction stage for ||v_k||
- *
- * All vectors v_i are stored column-major in the array 'device_krylov_basis', where each vector has an actual length 'v_k_size', but might be padded to have 'v_k_internal_size'
- */
-template <typename T>
-void pipelined_gmres_gram_schmidt_stage2(vector_base<T> & device_krylov_basis,
- vcl_size_t v_k_size,
- vcl_size_t v_k_internal_size,
- vcl_size_t k,
- vector_base<T> const & vi_in_vk_buffer,
- vector_base<T> & R_buffer,
- vcl_size_t krylov_dim,
- vector_base<T> & inner_prod_buffer,
- vcl_size_t buffer_chunk_size)
-{
- switch (viennacl::traits::handle(device_krylov_basis).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_gmres_gram_schmidt_stage2(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, R_buffer, krylov_dim, inner_prod_buffer, buffer_chunk_size);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_gmres_gram_schmidt_stage2(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, R_buffer, krylov_dim, inner_prod_buffer, buffer_chunk_size);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_gmres_gram_schmidt_stage2(device_krylov_basis, v_k_size, v_k_internal_size, k, vi_in_vk_buffer, R_buffer, krylov_dim, inner_prod_buffer, buffer_chunk_size);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-/** @brief Computes x += eta_0 r + sum_{i=1}^{k-1} eta_i v_{i-1} */
-template <typename T>
-void pipelined_gmres_update_result(vector_base<T> & result,
- vector_base<T> const & residual,
- vector_base<T> const & krylov_basis,
- vcl_size_t v_k_size,
- vcl_size_t v_k_internal_size,
- vector_base<T> const & coefficients,
- vcl_size_t k)
-{
- switch (viennacl::traits::handle(result).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_gmres_update_result(result, residual, krylov_basis, v_k_size, v_k_internal_size, coefficients, k);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_gmres_update_result(result, residual, krylov_basis, v_k_size, v_k_internal_size, coefficients, k);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_gmres_update_result(result, residual, krylov_basis, v_k_size, v_k_internal_size, coefficients, k);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-/** @brief Performs a joint vector update operation needed for an efficient pipelined GMRES algorithm.
- *
- * This routines computes for a matrix A and vectors 'p' and 'Ap':
- * Ap = prod(A, p);
- * and computes the two reduction stages for computing inner_prod(p,Ap), inner_prod(Ap,Ap)
- */
-template <typename MatrixType, typename T>
-void pipelined_gmres_prod(MatrixType const & A,
- vector_base<T> const & p,
- vector_base<T> & Ap,
- vector_base<T> & inner_prod_buffer)
-{
- switch (viennacl::traits::handle(p).get_active_handle_id())
- {
- case viennacl::MAIN_MEMORY:
- viennacl::linalg::host_based::pipelined_gmres_prod(A, p, Ap, inner_prod_buffer);
- break;
-#ifdef VIENNACL_WITH_OPENCL
- case viennacl::OPENCL_MEMORY:
- viennacl::linalg::opencl::pipelined_gmres_prod(A, p, Ap, inner_prod_buffer);
- break;
-#endif
-#ifdef VIENNACL_WITH_CUDA
- case viennacl::CUDA_MEMORY:
- viennacl::linalg::cuda::pipelined_gmres_prod(A, p, Ap, inner_prod_buffer);
- break;
-#endif
- case viennacl::MEMORY_NOT_INITIALIZED:
- throw memory_exception("not initialised!");
- default:
- throw memory_exception("not implemented");
- }
-}
-
-
-} //namespace linalg
-} //namespace viennacl
-
-
-#endif
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/jacobi_precond.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/jacobi_precond.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/jacobi_precond.hpp
deleted file mode 100644
index 0b16964..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/jacobi_precond.hpp
+++ /dev/null
@@ -1,141 +0,0 @@
-#ifndef VIENNACL_LINALG_JACOBI_PRECOND_HPP_
-#define VIENNACL_LINALG_JACOBI_PRECOND_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/jacobi_precond.hpp
- @brief Implementation of a simple Jacobi preconditioner
-*/
-
-#include <vector>
-#include <cmath>
-#include "viennacl/forwards.h"
-#include "viennacl/vector.hpp"
-#include "viennacl/compressed_matrix.hpp"
-#include "viennacl/tools/tools.hpp"
-#include "viennacl/linalg/sparse_matrix_operations.hpp"
-#include "viennacl/linalg/row_scaling.hpp"
-
-#include <map>
-
-namespace viennacl
-{
-namespace linalg
-{
-
-/** @brief A tag for a jacobi preconditioner
-*/
-class jacobi_tag {};
-
-
-/** @brief Jacobi preconditioner class, can be supplied to solve()-routines. Generic version for non-ViennaCL matrices.
-*/
-template<typename MatrixT,
- bool is_viennacl = detail::row_scaling_for_viennacl<MatrixT>::value >
-class jacobi_precond
-{
- typedef typename MatrixT::value_type NumericType;
-
- public:
- jacobi_precond(MatrixT const & mat, jacobi_tag const &) : diag_A_(viennacl::traits::size1(mat))
- {
- init(mat);
- }
-
- void init(MatrixT const & mat)
- {
- diag_A_.resize(viennacl::traits::size1(mat)); //resize without preserving values
-
- for (typename MatrixT::const_iterator1 row_it = mat.begin1();
- row_it != mat.end1();
- ++row_it)
- {
- bool diag_found = false;
- for (typename MatrixT::const_iterator2 col_it = row_it.begin();
- col_it != row_it.end();
- ++col_it)
- {
- if (col_it.index1() == col_it.index2())
- {
- diag_A_[col_it.index1()] = *col_it;
- diag_found = true;
- }
- }
- if (!diag_found)
- throw zero_on_diagonal_exception("ViennaCL: Zero in diagonal encountered while setting up Jacobi preconditioner!");
- }
- }
-
-
- /** @brief Apply to res = b - Ax, i.e. jacobi applied vec (right hand side), */
- template<typename VectorT>
- void apply(VectorT & vec) const
- {
- assert(viennacl::traits::size(diag_A_) == viennacl::traits::size(vec) && bool("Size mismatch"));
- for (vcl_size_t i=0; i<diag_A_.size(); ++i)
- vec[i] /= diag_A_[i];
- }
-
- private:
- std::vector<NumericType> diag_A_;
-};
-
-
-/** @brief Jacobi preconditioner class, can be supplied to solve()-routines.
-*
-* Specialization for compressed_matrix
-*/
-template<typename MatrixT>
-class jacobi_precond<MatrixT, true>
-{
- typedef typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type NumericType;
-
- public:
- jacobi_precond(MatrixT const & mat, jacobi_tag const &) : diag_A_(mat.size1(), viennacl::traits::context(mat))
- {
- init(mat);
- }
-
-
- void init(MatrixT const & mat)
- {
- detail::row_info(mat, diag_A_, detail::SPARSE_ROW_DIAGONAL);
- }
-
-
- template<unsigned int AlignmentV>
- void apply(viennacl::vector<NumericType, AlignmentV> & vec) const
- {
- assert(viennacl::traits::size(diag_A_) == viennacl::traits::size(vec) && bool("Size mismatch"));
- vec = element_div(vec, diag_A_);
- }
-
- private:
- viennacl::vector<NumericType> diag_A_;
-};
-
-}
-}
-
-
-
-
-#endif
-
-
-
http://git-wip-us.apache.org/repos/asf/mahout/blob/7ae549fa/native-viennaCL/src/main/cpp/viennacl/linalg/lanczos.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/lanczos.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/lanczos.hpp
deleted file mode 100644
index ffac471..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/lanczos.hpp
+++ /dev/null
@@ -1,515 +0,0 @@
-#ifndef VIENNACL_LINALG_LANCZOS_HPP_
-#define VIENNACL_LINALG_LANCZOS_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/lanczos.hpp
-* @brief Generic interface for the Lanczos algorithm.
-*
-* Contributed by Guenther Mader and Astrid Rupp.
-*/
-
-#include <cmath>
-#include <vector>
-#include "viennacl/vector.hpp"
-#include "viennacl/compressed_matrix.hpp"
-#include "viennacl/linalg/prod.hpp"
-#include "viennacl/linalg/inner_prod.hpp"
-#include "viennacl/linalg/norm_2.hpp"
-#include "viennacl/io/matrix_market.hpp"
-#include "viennacl/linalg/bisect.hpp"
-#include "viennacl/tools/random.hpp"
-
-namespace viennacl
-{
-namespace linalg
-{
-
-/** @brief A tag for the lanczos algorithm.
-*/
-class lanczos_tag
-{
-public:
-
- enum
- {
- partial_reorthogonalization = 0,
- full_reorthogonalization,
- no_reorthogonalization
- };
-
- /** @brief The constructor
- *
- * @param factor Exponent of epsilon - tolerance for batches of Reorthogonalization
- * @param numeig Number of eigenvalues to be returned
- * @param met Method for Lanczos-Algorithm: 0 for partial Reorthogonalization, 1 for full Reorthogonalization and 2 for Lanczos without Reorthogonalization
- * @param krylov Maximum krylov-space size
- */
-
- lanczos_tag(double factor = 0.75,
- vcl_size_t numeig = 10,
- int met = 0,
- vcl_size_t krylov = 100) : factor_(factor), num_eigenvalues_(numeig), method_(met), krylov_size_(krylov) {}
-
- /** @brief Sets the number of eigenvalues */
- void num_eigenvalues(vcl_size_t numeig){ num_eigenvalues_ = numeig; }
-
- /** @brief Returns the number of eigenvalues */
- vcl_size_t num_eigenvalues() const { return num_eigenvalues_; }
-
- /** @brief Sets the exponent of epsilon. Values between 0.6 and 0.9 usually give best results. */
- void factor(double fct) { factor_ = fct; }
-
- /** @brief Returns the exponent */
- double factor() const { return factor_; }
-
- /** @brief Sets the size of the kylov space. Must be larger than number of eigenvalues to compute. */
- void krylov_size(vcl_size_t max) { krylov_size_ = max; }
-
- /** @brief Returns the size of the kylov space */
- vcl_size_t krylov_size() const { return krylov_size_; }
-
- /** @brief Sets the reorthogonalization method */
- void method(int met){ method_ = met; }
-
- /** @brief Returns the reorthogonalization method */
- int method() const { return method_; }
-
-
-private:
- double factor_;
- vcl_size_t num_eigenvalues_;
- int method_; // see enum defined above for possible values
- vcl_size_t krylov_size_;
-};
-
-
-namespace detail
-{
- /** @brief Inverse iteration for finding an eigenvector for an eigenvalue.
- *
- * beta[0] to be ignored for consistency.
- */
- template<typename NumericT>
- void inverse_iteration(std::vector<NumericT> const & alphas, std::vector<NumericT> const & betas,
- NumericT & eigenvalue, std::vector<NumericT> & eigenvector)
- {
- std::vector<NumericT> alpha_sweeped = alphas;
- for (vcl_size_t i=0; i<alpha_sweeped.size(); ++i)
- alpha_sweeped[i] -= eigenvalue;
- for (vcl_size_t row=1; row < alpha_sweeped.size(); ++row)
- alpha_sweeped[row] -= betas[row] * betas[row] / alpha_sweeped[row-1];
-
- // starting guess: ignore last equation
- eigenvector[alphas.size() - 1] = 1.0;
-
- for (vcl_size_t iter=0; iter<1; ++iter)
- {
- // solve first n-1 equations (A - \lambda I) y = -beta[n]
- eigenvector[alphas.size() - 1] /= alpha_sweeped[alphas.size() - 1];
- for (vcl_size_t row2=1; row2 < alphas.size(); ++row2)
- {
- vcl_size_t row = alphas.size() - row2 - 1;
- eigenvector[row] -= eigenvector[row+1] * betas[row+1];
- eigenvector[row] /= alpha_sweeped[row];
- }
-
- // normalize eigenvector:
- NumericT norm_vector = 0;
- for (vcl_size_t i=0; i<eigenvector.size(); ++i)
- norm_vector += eigenvector[i] * eigenvector[i];
- norm_vector = std::sqrt(norm_vector);
- for (vcl_size_t i=0; i<eigenvector.size(); ++i)
- eigenvector[i] /= norm_vector;
- }
-
- //eigenvalue = (alphas[0] * eigenvector[0] + betas[1] * eigenvector[1]) / eigenvector[0];
- }
-
- /**
- * @brief Implementation of the Lanczos PRO algorithm (partial reorthogonalization)
- *
- * @param A The system matrix
- * @param r Random start vector
- * @param eigenvectors_A Dense matrix holding the eigenvectors of A (one eigenvector per column)
- * @param size Size of krylov-space
- * @param tag Lanczos_tag with several options for the algorithm
- * @param compute_eigenvectors Boolean flag. If true, eigenvectors are computed. Otherwise the routine returns after calculating eigenvalues.
- * @return Returns the eigenvalues (number of eigenvalues equals size of krylov-space)
- */
-
- template<typename MatrixT, typename DenseMatrixT, typename NumericT>
- std::vector<NumericT>
- lanczosPRO (MatrixT const& A, vector_base<NumericT> & r, DenseMatrixT & eigenvectors_A, vcl_size_t size, lanczos_tag const & tag, bool compute_eigenvectors)
- {
- // generation of some random numbers, used for lanczos PRO algorithm
- viennacl::tools::normal_random_numbers<NumericT> get_N;
-
- std::vector<vcl_size_t> l_bound(size/2), u_bound(size/2);
- vcl_size_t n = r.size();
- std::vector<NumericT> w(size), w_old(size); //w_k, w_{k-1}
-
- NumericT inner_rt;
- std::vector<NumericT> alphas, betas;
- viennacl::matrix<NumericT, viennacl::column_major> Q(n, size); //column-major matrix holding the Krylov basis vectors
-
- bool second_step = false;
- NumericT eps = std::numeric_limits<NumericT>::epsilon();
- NumericT squ_eps = std::sqrt(eps);
- NumericT eta = std::exp(std::log(eps) * tag.factor());
-
- NumericT beta = viennacl::linalg::norm_2(r);
-
- r /= beta;
-
- viennacl::vector_base<NumericT> q_0(Q.handle(), Q.size1(), 0, 1);
- q_0 = r;
-
- viennacl::vector<NumericT> u = viennacl::linalg::prod(A, r);
- NumericT alpha = viennacl::linalg::inner_prod(u, r);
- alphas.push_back(alpha);
- w[0] = 1;
- betas.push_back(beta);
-
- vcl_size_t batches = 0;
- for (vcl_size_t i = 1; i < size; i++) // Main loop for setting up the Krylov space
- {
- viennacl::vector_base<NumericT> q_iminus1(Q.handle(), Q.size1(), (i-1) * Q.internal_size1(), 1);
- r = u - alpha * q_iminus1;
- beta = viennacl::linalg::norm_2(r);
-
- betas.push_back(beta);
- r = r / beta;
-
- //
- // Update recurrence relation for estimating orthogonality loss
- //
- w_old = w;
- w[0] = (betas[1] * w_old[1] + (alphas[0] - alpha) * w_old[0] - betas[i - 1] * w_old[0]) / beta + eps * 0.3 * get_N() * (betas[1] + beta);
- for (vcl_size_t j = 1; j < i - 1; j++)
- w[j] = (betas[j + 1] * w_old[j + 1] + (alphas[j] - alpha) * w_old[j] + betas[j] * w_old[j - 1] - betas[i - 1] * w_old[j]) / beta + eps * 0.3 * get_N() * (betas[j + 1] + beta);
- w[i-1] = 0.6 * eps * NumericT(n) * get_N() * betas[1] / beta;
-
- //
- // Check whether there has been a need for reorthogonalization detected in the previous iteration.
- // If so, run the reorthogonalization for each batch
- //
- if (second_step)
- {
- for (vcl_size_t j = 0; j < batches; j++)
- {
- for (vcl_size_t k = l_bound[j] + 1; k < u_bound[j] - 1; k++)
- {
- viennacl::vector_base<NumericT> q_k(Q.handle(), Q.size1(), k * Q.internal_size1(), 1);
- inner_rt = viennacl::linalg::inner_prod(r, q_k);
- r = r - inner_rt * q_k;
- w[k] = 1.5 * eps * get_N();
- }
- }
- NumericT temp = viennacl::linalg::norm_2(r);
- r = r / temp;
- beta = beta * temp;
- second_step = false;
- }
- batches = 0;
-
- //
- // Check for semiorthogonality
- //
- for (vcl_size_t j = 0; j < i; j++)
- {
- if (std::fabs(w[j]) >= squ_eps) // tentative loss of orthonormality, hence reorthonomalize
- {
- viennacl::vector_base<NumericT> q_j(Q.handle(), Q.size1(), j * Q.internal_size1(), 1);
- inner_rt = viennacl::linalg::inner_prod(r, q_j);
- r = r - inner_rt * q_j;
- w[j] = 1.5 * eps * get_N();
- vcl_size_t k = j - 1;
-
- // orthogonalization with respect to earlier basis vectors
- while (std::fabs(w[k]) > eta)
- {
- viennacl::vector_base<NumericT> q_k(Q.handle(), Q.size1(), k * Q.internal_size1(), 1);
- inner_rt = viennacl::linalg::inner_prod(r, q_k);
- r = r - inner_rt * q_k;
- w[k] = 1.5 * eps * get_N();
- if (k == 0) break;
- k--;
- }
- l_bound[batches] = k;
-
- // orthogonalization with respect to later basis vectors
- k = j + 1;
- while (k < i && std::fabs(w[k]) > eta)
- {
- viennacl::vector_base<NumericT> q_k(Q.handle(), Q.size1(), k * Q.internal_size1(), 1);
- inner_rt = viennacl::linalg::inner_prod(r, q_k);
- r = r - inner_rt * q_k;
- w[k] = 1.5 * eps * get_N();
- k++;
- }
- u_bound[batches] = k - 1;
- batches++;
-
- j = k-1; // go to end of current batch
- }
- }
-
- //
- // Normalize basis vector and reorthogonalize as needed
- //
- if (batches > 0)
- {
- NumericT temp = viennacl::linalg::norm_2(r);
- r = r / temp;
- beta = beta * temp;
- second_step = true;
- }
-
- // store Krylov vector in Q:
- viennacl::vector_base<NumericT> q_i(Q.handle(), Q.size1(), i * Q.internal_size1(), 1);
- q_i = r;
-
- //
- // determine and store alpha = <r, u> with u = A q_i - beta q_{i-1}
- //
- u = viennacl::linalg::prod(A, r);
- u += (-beta) * q_iminus1;
- alpha = viennacl::linalg::inner_prod(u, r);
- alphas.push_back(alpha);
- }
-
- //
- // Step 2: Compute eigenvalues of tridiagonal matrix obtained during Lanczos iterations:
- //
- std::vector<NumericT> eigenvalues = bisect(alphas, betas);
-
- //
- // Step 3: Compute eigenvectors via inverse iteration. Does not update eigenvalues, so only approximate by nature.
- //
- if (compute_eigenvectors)
- {
- std::vector<NumericT> eigenvector_tridiag(alphas.size());
- for (std::size_t i=0; i < tag.num_eigenvalues(); ++i)
- {
- // compute eigenvector of tridiagonal matrix via inverse:
- inverse_iteration(alphas, betas, eigenvalues[eigenvalues.size() - i - 1], eigenvector_tridiag);
-
- // eigenvector w of full matrix A. Given as w = Q * u, where u is the eigenvector of the tridiagonal matrix
- viennacl::vector<NumericT> eigenvector_u(eigenvector_tridiag.size());
- viennacl::copy(eigenvector_tridiag, eigenvector_u);
-
- viennacl::vector_base<NumericT> eigenvector_A(eigenvectors_A.handle(),
- eigenvectors_A.size1(),
- eigenvectors_A.row_major() ? i : i * eigenvectors_A.internal_size1(),
- eigenvectors_A.row_major() ? eigenvectors_A.internal_size2() : 1);
- eigenvector_A = viennacl::linalg::prod(project(Q,
- range(0, Q.size1()),
- range(0, eigenvector_u.size())),
- eigenvector_u);
- }
- }
-
- return eigenvalues;
- }
-
-
- /**
- * @brief Implementation of the Lanczos FRO algorithm
- *
- * @param A The system matrix
- * @param r Random start vector
- * @param eigenvectors_A A dense matrix in which the eigenvectors of A will be stored. Both row- and column-major matrices are supported.
- * @param krylov_dim Size of krylov-space
- * @param tag The Lanczos tag holding tolerances, etc.
- * @param compute_eigenvectors Boolean flag. If true, eigenvectors are computed. Otherwise the routine returns after calculating eigenvalues.
- * @return Returns the eigenvalues (number of eigenvalues equals size of krylov-space)
- */
- template< typename MatrixT, typename DenseMatrixT, typename NumericT>
- std::vector<NumericT>
- lanczos(MatrixT const& A, vector_base<NumericT> & r, DenseMatrixT & eigenvectors_A, vcl_size_t krylov_dim, lanczos_tag const & tag, bool compute_eigenvectors)
- {
- std::vector<NumericT> alphas, betas;
- viennacl::vector<NumericT> Aq(r.size());
- viennacl::matrix<NumericT, viennacl::column_major> Q(r.size(), krylov_dim + 1); // Krylov basis (each Krylov vector is one column)
-
- NumericT norm_r = norm_2(r);
- NumericT beta = norm_r;
- r /= norm_r;
-
- // first Krylov vector:
- viennacl::vector_base<NumericT> q0(Q.handle(), Q.size1(), 0, 1);
- q0 = r;
-
- //
- // Step 1: Run Lanczos' method to obtain tridiagonal matrix
- //
- for (vcl_size_t i = 0; i < krylov_dim; i++)
- {
- betas.push_back(beta);
- // last available vector from Krylov basis:
- viennacl::vector_base<NumericT> q_i(Q.handle(), Q.size1(), i * Q.internal_size1(), 1);
-
- // Lanczos algorithm:
- // - Compute A * q:
- Aq = viennacl::linalg::prod(A, q_i);
-
- // - Form Aq <- Aq - <Aq, q_i> * q_i - beta * q_{i-1}, where beta is ||q_i|| before normalization in previous iteration
- NumericT alpha = viennacl::linalg::inner_prod(Aq, q_i);
- Aq -= alpha * q_i;
-
- if (i > 0)
- {
- viennacl::vector_base<NumericT> q_iminus1(Q.handle(), Q.size1(), (i-1) * Q.internal_size1(), 1);
- Aq -= beta * q_iminus1;
-
- // Extra measures for improved numerical stability?
- if (tag.method() == lanczos_tag::full_reorthogonalization)
- {
- // Gram-Schmidt (re-)orthogonalization:
- // TODO: Reuse fast (pipelined) routines from GMRES or GEMV
- for (vcl_size_t j = 0; j < i; j++)
- {
- viennacl::vector_base<NumericT> q_j(Q.handle(), Q.size1(), j * Q.internal_size1(), 1);
- NumericT inner_rq = viennacl::linalg::inner_prod(Aq, q_j);
- Aq -= inner_rq * q_j;
- }
- }
- }
-
- // normalize Aq and add to Krylov basis at column i+1 in Q:
- beta = viennacl::linalg::norm_2(Aq);
- viennacl::vector_base<NumericT> q_iplus1(Q.handle(), Q.size1(), (i+1) * Q.internal_size1(), 1);
- q_iplus1 = Aq / beta;
-
- alphas.push_back(alpha);
- }
-
- //
- // Step 2: Compute eigenvalues of tridiagonal matrix obtained during Lanczos iterations:
- //
- std::vector<NumericT> eigenvalues = bisect(alphas, betas);
-
- //
- // Step 3: Compute eigenvectors via inverse iteration. Does not update eigenvalues, so only approximate by nature.
- //
- if (compute_eigenvectors)
- {
- std::vector<NumericT> eigenvector_tridiag(alphas.size());
- for (std::size_t i=0; i < tag.num_eigenvalues(); ++i)
- {
- // compute eigenvector of tridiagonal matrix via inverse:
- inverse_iteration(alphas, betas, eigenvalues[eigenvalues.size() - i - 1], eigenvector_tridiag);
-
- // eigenvector w of full matrix A. Given as w = Q * u, where u is the eigenvector of the tridiagonal matrix
- viennacl::vector<NumericT> eigenvector_u(eigenvector_tridiag.size());
- viennacl::copy(eigenvector_tridiag, eigenvector_u);
-
- viennacl::vector_base<NumericT> eigenvector_A(eigenvectors_A.handle(),
- eigenvectors_A.size1(),
- eigenvectors_A.row_major() ? i : i * eigenvectors_A.internal_size1(),
- eigenvectors_A.row_major() ? eigenvectors_A.internal_size2() : 1);
- eigenvector_A = viennacl::linalg::prod(project(Q,
- range(0, Q.size1()),
- range(0, eigenvector_u.size())),
- eigenvector_u);
- }
- }
-
- return eigenvalues;
- }
-
-} // end namespace detail
-
-/**
-* @brief Implementation of the calculation of eigenvalues using lanczos (with and without reorthogonalization).
-*
-* Implementation of Lanczos with partial reorthogonalization is implemented separately.
-*
-* @param matrix The system matrix
-* @param eigenvectors_A A dense matrix in which the eigenvectors of A will be stored. Both row- and column-major matrices are supported.
-* @param tag Tag with several options for the lanczos algorithm
-* @param compute_eigenvectors Boolean flag. If true, eigenvectors are computed. Otherwise the routine returns after calculating eigenvalues.
-* @return Returns the n largest eigenvalues (n defined in the lanczos_tag)
-*/
-template<typename MatrixT, typename DenseMatrixT>
-std::vector< typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type >
-eig(MatrixT const & matrix, DenseMatrixT & eigenvectors_A, lanczos_tag const & tag, bool compute_eigenvectors = true)
-{
- typedef typename viennacl::result_of::value_type<MatrixT>::type NumericType;
- typedef typename viennacl::result_of::cpu_value_type<NumericType>::type CPU_NumericType;
- typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type VectorT;
-
- viennacl::tools::uniform_random_numbers<CPU_NumericType> random_gen;
-
- std::vector<CPU_NumericType> eigenvalues;
- vcl_size_t matrix_size = matrix.size1();
- VectorT r(matrix_size);
- std::vector<CPU_NumericType> s(matrix_size);
-
- for (vcl_size_t i=0; i<s.size(); ++i)
- s[i] = CPU_NumericType(0.5) + random_gen();
-
- detail::copy_vec_to_vec(s,r);
-
- vcl_size_t size_krylov = (matrix_size < tag.krylov_size()) ? matrix_size
- : tag.krylov_size();
-
- switch (tag.method())
- {
- case lanczos_tag::partial_reorthogonalization:
- eigenvalues = detail::lanczosPRO(matrix, r, eigenvectors_A, size_krylov, tag, compute_eigenvectors);
- break;
- case lanczos_tag::full_reorthogonalization:
- case lanczos_tag::no_reorthogonalization:
- eigenvalues = detail::lanczos(matrix, r, eigenvectors_A, size_krylov, tag, compute_eigenvectors);
- break;
- }
-
- std::vector<CPU_NumericType> largest_eigenvalues;
-
- for (vcl_size_t i = 1; i<=tag.num_eigenvalues(); i++)
- largest_eigenvalues.push_back(eigenvalues[size_krylov-i]);
-
-
- return largest_eigenvalues;
-}
-
-
-/**
-* @brief Implementation of the calculation of eigenvalues using lanczos (with and without reorthogonalization).
-*
-* Implementation of Lanczos with partial reorthogonalization is implemented separately.
-*
-* @param matrix The system matrix
-* @param tag Tag with several options for the lanczos algorithm
-* @return Returns the n largest eigenvalues (n defined in the lanczos_tag)
-*/
-template<typename MatrixT>
-std::vector< typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type >
-eig(MatrixT const & matrix, lanczos_tag const & tag)
-{
- typedef typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type NumericType;
-
- viennacl::matrix<NumericType> eigenvectors(matrix.size1(), tag.num_eigenvalues());
- return eig(matrix, eigenvectors, tag, false);
-}
-
-} // 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/lu.hpp
----------------------------------------------------------------------
diff --git a/native-viennaCL/src/main/cpp/viennacl/linalg/lu.hpp b/native-viennaCL/src/main/cpp/viennacl/linalg/lu.hpp
deleted file mode 100644
index 0bdd037..0000000
--- a/native-viennaCL/src/main/cpp/viennacl/linalg/lu.hpp
+++ /dev/null
@@ -1,227 +0,0 @@
-#ifndef VIENNACL_LINALG_LU_HPP
-#define VIENNACL_LINALG_LU_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/lu.hpp
- @brief Implementations of LU factorization for row-major and column-major dense matrices.
-*/
-
-#include <algorithm> //for std::min
-
-#include "viennacl/matrix.hpp"
-#include "viennacl/matrix_proxy.hpp"
-
-#include "viennacl/linalg/prod.hpp"
-#include "viennacl/linalg/direct_solve.hpp"
-
-namespace viennacl
-{
-namespace linalg
-{
-/** @brief LU factorization of a row-major dense matrix.
-*
-* @param A The system matrix, where the LU matrices are directly written to. The implicit unit diagonal of L is not written.
-*/
-template<typename NumericT>
-void lu_factorize(matrix<NumericT, viennacl::row_major> & A)
-{
- typedef matrix<NumericT, viennacl::row_major> MatrixType;
-
- vcl_size_t max_block_size = 32;
- vcl_size_t num_blocks = (A.size2() - 1) / max_block_size + 1;
- std::vector<NumericT> temp_buffer(A.internal_size2() * max_block_size);
-
- // Iterate over panels
- for (vcl_size_t panel_id = 0; panel_id < num_blocks; ++panel_id)
- {
- vcl_size_t row_start = panel_id * max_block_size;
- vcl_size_t current_block_size = std::min<vcl_size_t>(A.size1() - row_start, max_block_size);
-
- viennacl::range block_range(row_start, row_start + current_block_size);
- viennacl::range remainder_range(row_start + current_block_size, A.size1());
-
- //
- // Perform LU factorization on panel:
- //
-
-
- // Read from matrix to buffer:
- viennacl::backend::memory_read(A.handle(),
- sizeof(NumericT) * row_start * A.internal_size2(),
- sizeof(NumericT) * current_block_size * A.internal_size2(),
- &(temp_buffer[0]));
-
- // Factorize (kij-version):
- for (vcl_size_t k=0; k < current_block_size - 1; ++k)
- {
- for (vcl_size_t i=k+1; i < current_block_size; ++i)
- {
- temp_buffer[row_start + i * A.internal_size2() + k] /= temp_buffer[row_start + k * A.internal_size2() + k]; // write l_ik
-
- NumericT l_ik = temp_buffer[row_start + i * A.internal_size2() + k];
-
- for (vcl_size_t j = row_start + k + 1; j < A.size1(); ++j)
- temp_buffer[i * A.internal_size2() + j] -= l_ik * temp_buffer[k * A.internal_size2() + j]; // l_ik * a_kj
- }
- }
-
- // Write back:
- viennacl::backend::memory_write(A.handle(),
- sizeof(NumericT) * row_start * A.internal_size2(),
- sizeof(NumericT) * current_block_size * A.internal_size2(),
- &(temp_buffer[0]));
-
- if (remainder_range.size() > 0)
- {
- //
- // Compute L_12 = [ (U_11)^{T}^{-1} A_{21}^T ]^T
- //
- viennacl::matrix_range<MatrixType> U_11(A, block_range, block_range);
- viennacl::matrix_range<MatrixType> A_21(A, remainder_range, block_range);
- viennacl::linalg::inplace_solve(trans(U_11), trans(A_21), viennacl::linalg::lower_tag());
-
- //
- // Update remainder of A
- //
- viennacl::matrix_range<MatrixType> L_21(A, remainder_range, block_range);
- viennacl::matrix_range<MatrixType> U_12(A, block_range, remainder_range);
- viennacl::matrix_range<MatrixType> A_22(A, remainder_range, remainder_range);
-
- A_22 -= viennacl::linalg::prod(L_21, U_12);
- }
- }
-
-}
-
-
-/** @brief LU factorization of a column-major dense matrix.
-*
-* @param A The system matrix, where the LU matrices are directly written to. The implicit unit diagonal of L is not written.
-*/
-template<typename NumericT>
-void lu_factorize(matrix<NumericT, viennacl::column_major> & A)
-{
- typedef matrix<NumericT, viennacl::column_major> MatrixType;
-
- vcl_size_t max_block_size = 32;
- vcl_size_t num_blocks = (A.size1() - 1) / max_block_size + 1;
- std::vector<NumericT> temp_buffer(A.internal_size1() * max_block_size);
-
- // Iterate over panels
- for (vcl_size_t panel_id = 0; panel_id < num_blocks; ++panel_id)
- {
- vcl_size_t col_start = panel_id * max_block_size;
- vcl_size_t current_block_size = std::min<vcl_size_t>(A.size1() - col_start, max_block_size);
-
- viennacl::range block_range(col_start, col_start + current_block_size);
- viennacl::range remainder_range(col_start + current_block_size, A.size1());
-
- //
- // Perform LU factorization on panel:
- //
-
-
- // Read from matrix to buffer:
- viennacl::backend::memory_read(A.handle(),
- sizeof(NumericT) * col_start * A.internal_size1(),
- sizeof(NumericT) * current_block_size * A.internal_size1(),
- &(temp_buffer[0]));
-
- // Factorize (kji-version):
- for (vcl_size_t k=0; k < current_block_size; ++k)
- {
- NumericT a_kk = temp_buffer[col_start + k + k * A.internal_size1()];
- for (vcl_size_t i=col_start+k+1; i < A.size1(); ++i)
- temp_buffer[i + k * A.internal_size1()] /= a_kk; // write l_ik
-
- for (vcl_size_t j=k+1; j < current_block_size; ++j)
- {
- NumericT a_kj = temp_buffer[col_start + k + j * A.internal_size1()];
- for (vcl_size_t i=col_start+k+1; i < A.size1(); ++i)
- temp_buffer[i + j * A.internal_size1()] -= temp_buffer[i + k * A.internal_size1()] * a_kj; // l_ik * a_kj
- }
- }
-
- // Write back:
- viennacl::backend::memory_write(A.handle(),
- sizeof(NumericT) * col_start * A.internal_size1(),
- sizeof(NumericT) * current_block_size * A.internal_size1(),
- &(temp_buffer[0]));
-
- if (remainder_range.size() > 0)
- {
- //
- // Compute U_12:
- //
- viennacl::matrix_range<MatrixType> L_11(A, block_range, block_range);
- viennacl::matrix_range<MatrixType> A_12(A, block_range, remainder_range);
- viennacl::linalg::inplace_solve(L_11, A_12, viennacl::linalg::unit_lower_tag());
-
- //
- // Update remainder of A
- //
- viennacl::matrix_range<MatrixType> L_21(A, remainder_range, block_range);
- viennacl::matrix_range<MatrixType> U_12(A, block_range, remainder_range);
- viennacl::matrix_range<MatrixType> A_22(A, remainder_range, remainder_range);
-
- A_22 -= viennacl::linalg::prod(L_21, U_12);
- }
-
- }
-
-}
-
-
-//
-// Convenience layer:
-//
-
-/** @brief LU substitution for the system LU = rhs.
-*
-* @param A The system matrix, where the LU matrices are directly written to. The implicit unit diagonal of L is not written.
-* @param B The matrix of load vectors, where the solution is directly written to
-*/
-template<typename NumericT, typename F1, typename F2, unsigned int AlignmentV1, unsigned int AlignmentV2>
-void lu_substitute(matrix<NumericT, F1, AlignmentV1> const & A,
- matrix<NumericT, F2, AlignmentV2> & B)
-{
- assert(A.size1() == A.size2() && bool("Matrix must be square"));
- assert(A.size1() == B.size1() && bool("Matrix must be square"));
- inplace_solve(A, B, unit_lower_tag());
- inplace_solve(A, B, upper_tag());
-}
-
-/** @brief LU substitution for the system LU = rhs.
-*
-* @param A The system matrix, where the LU matrices are directly written to. The implicit unit diagonal of L is not written.
-* @param vec The load vector, where the solution is directly written to
-*/
-template<typename NumericT, typename F, unsigned int MatAlignmentV, unsigned int VecAlignmentV>
-void lu_substitute(matrix<NumericT, F, MatAlignmentV> const & A,
- vector<NumericT, VecAlignmentV> & vec)
-{
- assert(A.size1() == A.size2() && bool("Matrix must be square"));
- inplace_solve(A, vec, unit_lower_tag());
- inplace_solve(A, vec, upper_tag());
-}
-
-}
-}
-
-#endif