All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ilu_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_OPENCL_ILU_OPERATIONS_HPP_
#define VIENNACL_LINALG_OPENCL_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
============================================================================= */
#include <cmath>
#include <algorithm> //for std::max and std::min
#include "viennacl/forwards.h"
#include "viennacl/scalar.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/opencl/common.hpp"
#include "viennacl/linalg/opencl/kernels/ilu.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/vector_operations.hpp"
namespace viennacl
{
namespace linalg
{
namespace opencl
{
template<typename NumericT>
void extract_L(compressed_matrix<NumericT> const & A,
compressed_matrix<NumericT> & L)
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
//
// Step 1: Count elements in L:
//
viennacl::ocl::kernel & k1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "extract_L_1");
viennacl::ocl::enqueue(k1(A.handle1().opencl_handle(), A.handle2().opencl_handle(), cl_uint(A.size1()),
L.handle1().opencl_handle())
);
//
// Step 2: Exclusive scan on row_buffers:
//
viennacl::vector_base<unsigned int> wrapped_L_row_buffer(L.handle1(), A.size1() + 1, 0, 1);
viennacl::linalg::exclusive_scan(wrapped_L_row_buffer, wrapped_L_row_buffer);
L.reserve(wrapped_L_row_buffer[L.size1()], false);
//
// Step 3: Write entries
//
viennacl::ocl::kernel & k2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "extract_L_2");
viennacl::ocl::enqueue(k2(A.handle1().opencl_handle(), A.handle2().opencl_handle(), A.handle().opencl_handle(), cl_uint(A.size1()),
L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle())
);
L.generate_row_block_information();
} // extract_LU
template<typename NumericT>
void icc_scale(compressed_matrix<NumericT> const & A,
compressed_matrix<NumericT> & L)
{
viennacl::vector<NumericT> D(A.size1(), viennacl::traits::context(A));
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
// fill D:
viennacl::ocl::kernel & k1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_scale_kernel_1");
viennacl::ocl::enqueue(k1(A.handle1().opencl_handle(), A.handle2().opencl_handle(), A.handle().opencl_handle(), cl_uint(A.size1()), D) );
// scale L:
viennacl::ocl::kernel & k2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_scale_kernel_2");
viennacl::ocl::enqueue(k2(L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle(), cl_uint(A.size1()), D) );
}
template<typename NumericT>
void icc_chow_patel_sweep(compressed_matrix<NumericT> & L,
vector<NumericT> const & aij_L)
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(L).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
viennacl::backend::mem_handle L_backup;
viennacl::backend::memory_create(L_backup, L.handle().raw_size(), viennacl::traits::context(L));
viennacl::backend::memory_copy(L.handle(), L_backup, 0, 0, L.handle().raw_size());
viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "icc_chow_patel_sweep_kernel");
viennacl::ocl::enqueue(k(L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle(), L_backup.opencl_handle(), cl_uint(L.size1()),
aij_L)
);
}
template<typename NumericT>
void extract_LU(compressed_matrix<NumericT> const & A,
compressed_matrix<NumericT> & L,
compressed_matrix<NumericT> & U)
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
//
// Step 1: Count elements in L and U:
//
viennacl::ocl::kernel & k1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "extract_LU_1");
viennacl::ocl::enqueue(k1(A.handle1().opencl_handle(), A.handle2().opencl_handle(), cl_uint(A.size1()),
L.handle1().opencl_handle(),
U.handle1().opencl_handle())
);
//
// Step 2: Exclusive scan on row_buffers:
//
viennacl::vector_base<unsigned int> wrapped_L_row_buffer(L.handle1(), A.size1() + 1, 0, 1);
viennacl::linalg::exclusive_scan(wrapped_L_row_buffer, wrapped_L_row_buffer);
L.reserve(wrapped_L_row_buffer[L.size1()], false);
viennacl::vector_base<unsigned int> wrapped_U_row_buffer(U.handle1(), A.size1() + 1, 0, 1);
viennacl::linalg::exclusive_scan(wrapped_U_row_buffer, wrapped_U_row_buffer);
U.reserve(wrapped_U_row_buffer[U.size1()], false);
//
// Step 3: Write entries
//
viennacl::ocl::kernel & k2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "extract_LU_2");
viennacl::ocl::enqueue(k2(A.handle1().opencl_handle(), A.handle2().opencl_handle(), A.handle().opencl_handle(), cl_uint(A.size1()),
L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle(),
U.handle1().opencl_handle(), U.handle2().opencl_handle(), U.handle().opencl_handle())
);
L.generate_row_block_information();
// Note: block information for U will be generated after transposition
} // extract_LU
template<typename NumericT>
void ilu_scale(compressed_matrix<NumericT> const & A,
compressed_matrix<NumericT> & L,
compressed_matrix<NumericT> & U)
{
viennacl::vector<NumericT> D(A.size1(), viennacl::traits::context(A));
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
// fill D:
viennacl::ocl::kernel & k1 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_scale_kernel_1");
viennacl::ocl::enqueue(k1(A.handle1().opencl_handle(), A.handle2().opencl_handle(), A.handle().opencl_handle(), cl_uint(A.size1()), D) );
// scale L:
viennacl::ocl::kernel & k2 = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_scale_kernel_2");
viennacl::ocl::enqueue(k2(L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle(), cl_uint(A.size1()), D) );
// scale U:
viennacl::ocl::enqueue(k2(U.handle1().opencl_handle(), U.handle2().opencl_handle(), U.handle().opencl_handle(), cl_uint(A.size1()), D) );
}
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)
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(L).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
viennacl::backend::mem_handle L_backup;
viennacl::backend::memory_create(L_backup, L.handle().raw_size(), viennacl::traits::context(L));
viennacl::backend::memory_copy(L.handle(), L_backup, 0, 0, L.handle().raw_size());
viennacl::backend::mem_handle U_backup;
viennacl::backend::memory_create(U_backup, U_trans.handle().raw_size(), viennacl::traits::context(U_trans));
viennacl::backend::memory_copy(U_trans.handle(), U_backup, 0, 0, U_trans.handle().raw_size());
viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_chow_patel_sweep_kernel");
viennacl::ocl::enqueue(k(L.handle1().opencl_handle(), L.handle2().opencl_handle(), L.handle().opencl_handle(), L_backup.opencl_handle(), cl_uint(L.size1()),
aij_L,
U_trans.handle1().opencl_handle(), U_trans.handle2().opencl_handle(), U_trans.handle().opencl_handle(), U_backup.opencl_handle(),
aij_U_trans)
);
}
template<typename NumericT>
void ilu_form_neumann_matrix(compressed_matrix<NumericT> & R,
vector<NumericT> & diag_R)
{
viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(R).context());
viennacl::linalg::opencl::kernels::ilu<NumericT>::init(ctx);
viennacl::ocl::kernel & k = ctx.get_kernel(viennacl::linalg::opencl::kernels::ilu<NumericT>::program_name(), "ilu_form_neumann_matrix_kernel");
viennacl::ocl::enqueue(k(R.handle1().opencl_handle(), R.handle2().opencl_handle(), R.handle().opencl_handle(), cl_uint(R.size1()),
diag_R)
);
}
} //namespace opencl
} //namespace linalg
} //namespace viennacl
#endif
size.hpp
Generic size and resize functionality for different vector and matrix types.
vector_operations.hpp
Implementations of vector operations.
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
viennacl::ocl::kernel
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:58
start.hpp
Extracts the underlying OpenCL start index handle from a vector, a matrix, an expression etc...
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::ocl::context
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:55
viennacl::linalg::opencl::ilu_scale
void ilu_scale(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L, compressed_matrix< NumericT > &U)
Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L a...
Definition: ilu_operations.hpp:186
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
stride.hpp
Determines row and column increments for matrices and matrix proxies.
viennacl::linalg::opencl::ilu_chow_patel_sweep
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)
Performs one nonlinear relaxation step in the Chow-Patel-ILU using OpenCL (cf. Algorithm 2 in paper) ...
Definition: ilu_operations.hpp:213
viennacl::compressed_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:942
viennacl::compressed_matrix::handle1
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:936
common.hpp
Common implementations shared by OpenCL-based operations.
viennacl::compressed_matrix::generate_row_block_information
void generate_row_block_information()
Builds the row block information needed for fast sparse matrix-vector multiplications.
Definition: compressed_matrix.hpp:999
viennacl::linalg::opencl::kernels::ilu
Main kernel class for generating OpenCL kernels for incomplete LU factorization preconditioners.
Definition: ilu.hpp:453
viennacl::linalg::opencl::extract_L
void extract_L(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L)
Definition: ilu_operations.hpp:51
viennacl::ocl::context::get_kernel
viennacl::ocl::kernel & get_kernel(std::string const &program_name, std::string const &kernel_name)
Convenience function for retrieving the kernel of a program directly from the context.
Definition: context.hpp:605
ilu.hpp
OpenCL kernel file for nonnegative matrix factorization.
viennacl::linalg::opencl::kernels::ilu::init
static void init(viennacl::ocl::context &ctx)
Definition: ilu.hpp:460
viennacl::compressed_matrix::handle2
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:938
viennacl::vector_base< unsigned int >
viennacl::vector< NumericT >
predicate.hpp
All the predicates used within ViennaCL. Checks for expressions to be vectors, etc.
viennacl::compressed_matrix::reserve
void reserve(vcl_size_t new_nonzeros, bool preserve=true)
Allocate memory for the supplied number of nonzeros in the matrix. Old values are preserved...
Definition: compressed_matrix.hpp:794
viennacl::backend::memory_copy
void memory_copy(mem_handle const &src_buffer, mem_handle &dst_buffer, vcl_size_t src_offset, vcl_size_t dst_offset, vcl_size_t bytes_to_copy)
Copies 'bytes_to_copy' bytes from address 'src_buffer + src_offset' to memory starting at address 'ds...
Definition: memory.hpp:140
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::ocl::enqueue
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:50
viennacl::linalg::opencl::ilu_form_neumann_matrix
void ilu_form_neumann_matrix(compressed_matrix< NumericT > &R, vector< NumericT > &diag_R)
Definition: ilu_operations.hpp:243
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::backend::mem_handle::raw_size
vcl_size_t raw_size() const
Returns the number of bytes of the currently active buffer.
Definition: mem_handle.hpp:230
viennacl::compressed_matrix< NumericT >
viennacl::linalg::exclusive_scan
void exclusive_scan(vector_base< NumericT > &vec1, vector_base< NumericT > &vec2)
This function implements an exclusive scan.
Definition: vector_operations.hpp:1240
viennacl::backend::memory_create
void memory_create(mem_handle &handle, vcl_size_t size_in_bytes, viennacl::context const &ctx, const void *host_ptr=NULL)
Creates an array of the specified size. If the second argument is provided, the buffer is initialized...
Definition: memory.hpp:87
viennacl::linalg::opencl::icc_chow_patel_sweep
void icc_chow_patel_sweep(compressed_matrix< NumericT > &L, vector< NumericT > const &aij_L)
Performs one nonlinear relaxation step in the Chow-Patel-ILU using OpenCL (cf. Algorithm 2 in paper) ...
Definition: ilu_operations.hpp:116
viennacl::linalg::opencl::extract_LU
void extract_LU(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L, compressed_matrix< NumericT > &U)
Definition: ilu_operations.hpp:137
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::linalg::opencl::icc_scale
void icc_scale(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L)
Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L a...
Definition: ilu_operations.hpp:93
enable_if.hpp
Simple enable-if variant that uses the SFINAE pattern.