All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sparse_matrix_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_CUDA_SPARSE_MATRIX_OPERATIONS_HPP_
#define VIENNACL_LINALG_CUDA_SPARSE_MATRIX_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
============================================================================= */
#include "viennacl/forwards.h"
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/cuda/common.hpp"
#include "viennacl/linalg/cuda/vector_operations.hpp"
#include "viennacl/linalg/cuda/sparse_matrix_operations_solve.hpp"
//#ifdef VIENNACL_WITH_SPGEMM_RMERGE
#include "viennacl/linalg/cuda/spgemm_rmerge.hpp"
//#else
// #include "viennacl/linalg/cuda/spgemm.hpp"
//#endif
namespace viennacl
{
namespace linalg
{
namespace cuda
{
//
// Compressed matrix
//
namespace detail
{
template<typename NumericT>
__global__ void csr_row_info_extractor_kernel(
const unsigned int * row_indices,
const unsigned int * column_indices,
const NumericT * elements,
NumericT * result,
unsigned int size,
unsigned int option)
{
for (unsigned int row = blockDim.x * blockIdx.x + threadIdx.x;
row < size;
row += gridDim.x * blockDim.x)
{
NumericT value = 0;
unsigned int row_end = row_indices[row+1];
switch (option)
{
case 0: //inf-norm
for (unsigned int i = row_indices[row]; i < row_end; ++i)
value = max(value, fabs(elements[i]));
break;
case 1: //1-norm
for (unsigned int i = row_indices[row]; i < row_end; ++i)
value += fabs(elements[i]);
break;
case 2: //2-norm
for (unsigned int i = row_indices[row]; i < row_end; ++i)
value += elements[i] * elements[i];
value = sqrt(value);
break;
case 3: //diagonal entry
for (unsigned int i = row_indices[row]; i < row_end; ++i)
{
if (column_indices[i] == row)
{
value = elements[i];
break;
}
}
break;
default:
break;
}
result[row] = value;
}
}
template<typename NumericT, unsigned int AligmentV>
void row_info(compressed_matrix<NumericT, AligmentV> const & mat,
vector_base<NumericT> & vec,
viennacl::linalg::detail::row_info_types info_selector)
{
csr_row_info_extractor_kernel<<<128, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(info_selector)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_row_info_extractor_kernel");
}
struct spmv_pure
{
template<typename NumericT>
__device__ static void apply(NumericT & result, NumericT alpha, NumericT Ax, NumericT beta) { result = Ax; }
};
struct spmv_alpha_beta
{
template<typename NumericT>
__device__ static void apply(NumericT & result, NumericT alpha, NumericT Ax, NumericT beta) { result = alpha * Ax + ((beta != 0) ? beta * result : 0); }
};
} //namespace detail
template<unsigned int SubWarpSizeV, typename AlphaBetaHandlerT, typename NumericT>
__global__ void compressed_matrix_vec_mul_kernel(
const unsigned int * row_indices,
const unsigned int * column_indices,
const NumericT * elements,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
unsigned int size_result,
NumericT beta)
{
__shared__ NumericT shared_elements[512];
const unsigned int id_in_row = threadIdx.x % SubWarpSizeV;
const unsigned int block_increment = blockDim.x * ((size_result - 1) / (gridDim.x * blockDim.x) + 1);
const unsigned int block_start = blockIdx.x * block_increment;
const unsigned int block_stop = min(block_start + block_increment, size_result);
for (unsigned int row = block_start + threadIdx.x / SubWarpSizeV;
row < block_stop;
row += blockDim.x / SubWarpSizeV)
{
NumericT dot_prod = NumericT(0);
unsigned int row_end = row_indices[row+1];
for (unsigned int i = row_indices[row] + id_in_row; i < row_end; i += SubWarpSizeV)
dot_prod += elements[i] * x[column_indices[i] * inc_x + start_x];
shared_elements[threadIdx.x] = dot_prod;
if (1 < SubWarpSizeV) shared_elements[threadIdx.x] += shared_elements[threadIdx.x ^ 1];
if (2 < SubWarpSizeV) shared_elements[threadIdx.x] += shared_elements[threadIdx.x ^ 2];
if (4 < SubWarpSizeV) shared_elements[threadIdx.x] += shared_elements[threadIdx.x ^ 4];
if (8 < SubWarpSizeV) shared_elements[threadIdx.x] += shared_elements[threadIdx.x ^ 8];
if (16 < SubWarpSizeV) shared_elements[threadIdx.x] += shared_elements[threadIdx.x ^ 16];
if (id_in_row == 0)
AlphaBetaHandlerT::apply(result[row * inc_result + start_result], alpha, shared_elements[threadIdx.x], beta);
}
}
template<typename AlphaBetaHandlerT, typename NumericT>
__global__ void compressed_matrix_vec_mul_adaptive_kernel(
const unsigned int * row_indices,
const unsigned int * column_indices,
const unsigned int * row_blocks,
const NumericT * elements,
unsigned int num_blocks,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
unsigned int size_result,
NumericT beta)
{
__shared__ NumericT shared_elements[1024];
for (unsigned int block_id = blockIdx.x; block_id < num_blocks; block_id += gridDim.x)
{
unsigned int row_start = row_blocks[block_id];
unsigned int row_stop = row_blocks[block_id + 1];
unsigned int element_start = row_indices[row_start];
unsigned int element_stop = row_indices[row_stop];
unsigned int rows_to_process = row_stop - row_start;
if (rows_to_process > 1) // CSR stream with one thread per row
{
// load to shared buffer:
for (unsigned int i = element_start + threadIdx.x; i < element_stop; i += blockDim.x)
shared_elements[i - element_start] = elements[i] * x[column_indices[i] * inc_x + start_x];
__syncthreads();
// use one thread per row to sum:
for (unsigned int row = row_start + threadIdx.x; row < row_stop; row += blockDim.x)
{
NumericT dot_prod = 0;
unsigned int thread_row_start = row_indices[row] - element_start;
unsigned int thread_row_stop = row_indices[row + 1] - element_start;
for (unsigned int i = thread_row_start; i < thread_row_stop; ++i)
dot_prod += shared_elements[i];
AlphaBetaHandlerT::apply(result[row * inc_result + start_result], alpha, dot_prod, beta);
}
}
// TODO here: Consider CSR vector for two to four rows (cf. OpenCL implementation. Experience on Fermi suggests that this may not be necessary)
else // CSR vector for a single row
{
// load and sum to shared buffer:
shared_elements[threadIdx.x] = 0;
for (unsigned int i = element_start + threadIdx.x; i < element_stop; i += blockDim.x)
shared_elements[threadIdx.x] += elements[i] * x[column_indices[i] * inc_x + start_x];
// reduction to obtain final result
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
shared_elements[threadIdx.x] += shared_elements[threadIdx.x+stride];
}
if (threadIdx.x == 0)
AlphaBetaHandlerT::apply(result[row_start * inc_result + start_result], alpha, shared_elements[0], beta);
}
__syncthreads(); // avoid race conditions
}
}
template<class NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::compressed_matrix<NumericT, AlignmentV> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
static bool first = true;
static bool is_maxwell = false;
// check whether the CUDA device is from the Maxwell family.
// Only run once, because the query to the backend takes about the same time as a kernel launch (~15us), thus being too expensive to query each time.
//
// Note: This might result in non-optimal kernels being selected if multiple Maxwell- and non-Maxwell GPUs are available in the system and devices are switched at runtime.
// However, this situation is certainly rare, hence the the benefits of this singleton outweigh the disadvantages encountered in such a corner case.
if (first)
{
cudaDeviceProp prop;
int device_index = 0;
cudaError_t err_flag = cudaGetDevice(&device_index);
if (err_flag == cudaSuccess)
{
err_flag = cudaGetDeviceProperties(&prop, device_index);
if (err_flag == cudaSuccess && prop.major >= 5)
is_maxwell = true;
}
first = false;
}
if (is_maxwell && double(mat.nnz()) / double(mat.size1()) > 6.4) // less than 10% of threads expected to idle
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
compressed_matrix_vec_mul_kernel<8, detail::spmv_alpha_beta, NumericT><<<512, 256>>>( // experience on a GTX 750 Ti suggests that 8 is a substantially better choice here
viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
else
compressed_matrix_vec_mul_kernel<8, detail::spmv_pure, NumericT><<<512, 256>>>( // experience on a GTX 750 Ti suggests that 8 is a substantially better choice here
viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_vec_mul_kernel");
}
else if (!is_maxwell && double(mat.nnz()) / double(mat.size1()) > 12.0) // less than 25% of threads expected to idle
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
compressed_matrix_vec_mul_kernel<16, detail::spmv_alpha_beta, NumericT><<<512, 256>>>( // Fermi and Kepler prefer 16 threads per row (half-warp)
viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
else
compressed_matrix_vec_mul_kernel<16, detail::spmv_pure, NumericT><<<512, 256>>>( // Fermi and Kepler prefer 16 threads per row (half-warp)
viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_vec_mul_kernel");
}
else
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
compressed_matrix_vec_mul_adaptive_kernel<detail::spmv_alpha_beta><<<512, 256>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<NumericT>(mat.handle()),
static_cast<unsigned int>(mat.blocks1()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
else
compressed_matrix_vec_mul_adaptive_kernel<detail::spmv_pure><<<512, 256>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<NumericT>(mat.handle()),
static_cast<unsigned int>(mat.blocks1()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_vec_mul_adaptive_kernel");
}
}
template<typename LayoutT>
struct mat_mult_matrix_index
{
static __device__ unsigned int apply(unsigned int i, unsigned int j,
unsigned int row_start, unsigned int row_inc,
unsigned int col_start, unsigned int col_inc,
unsigned int internal_rows, unsigned int internal_cols)
{
return (row_start + i * row_inc) * internal_cols + col_start + j * col_inc;
}
};
template<>
struct mat_mult_matrix_index<viennacl::column_major>
{
static __device__ unsigned int apply(unsigned int i, unsigned int j,
unsigned int row_start, unsigned int row_inc,
unsigned int col_start, unsigned int col_inc,
unsigned int internal_rows, unsigned int internal_cols)
{
return (row_start + i * row_inc) + (col_start + j * col_inc) * internal_rows;
}
};
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void compressed_matrix_d_mat_mul_kernel(
const unsigned int * sp_mat_row_indices,
const unsigned int * sp_mat_col_indices,
const NumericT * sp_mat_elements,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
for (unsigned int row = blockIdx.x; row < result_row_size; row += gridDim.x)
{
unsigned int row_start = sp_mat_row_indices[row];
unsigned int row_end = sp_mat_row_indices[row+1];
for ( unsigned int col = threadIdx.x; col < result_col_size; col += blockDim.x)
{
NumericT r = 0;
for (unsigned int k = row_start; k < row_end; k++)
{
unsigned int j = sp_mat_col_indices[k];
NumericT x = sp_mat_elements[k];
NumericT y = d_mat[ DMatIndexT::apply(j, col,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols) ];
r += x * y;
}
result[ResultIndexT::apply(row, col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = r;
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::compressed_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_base<NumericT> & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.row_major() && result.row_major())
{
compressed_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_mat_mul_kernel");
}
else if (d_mat.row_major() && !result.row_major())
{
compressed_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_mat_mul_kernel");
}
else if (!d_mat.row_major() && result.row_major())
{
compressed_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_mat_mul_kernel");
}
else
{
compressed_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_mat_mul_kernel");
}
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void compressed_matrix_d_tr_mat_mul_kernel(
const unsigned int * sp_mat_row_indices,
const unsigned int * sp_mat_col_indices,
const NumericT * sp_mat_elements,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
for (unsigned int row = blockIdx.x; row < result_row_size; row += gridDim.x)
{
unsigned int row_start = sp_mat_row_indices[row];
unsigned int row_end = sp_mat_row_indices[row+1];
for ( unsigned int col = threadIdx.x; col < result_col_size; col += blockDim.x)
{
NumericT r = 0;
for (unsigned int k = row_start; k < row_end; k++)
{
unsigned int j = sp_mat_col_indices[k];
NumericT x = sp_mat_elements[k];
NumericT y = d_mat[ DMatIndexT::apply(col, j,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols) ];
r += x * y;
}
result [ ResultIndexT::apply(row, col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = r;
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::compressed_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
const viennacl::matrix_base<NumericT>,
viennacl::op_trans > & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.lhs().row_major() && result.row_major())
{
compressed_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_tr_mat_mul_kernel");
}
else if (d_mat.lhs().row_major() && !result.row_major())
{
compressed_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_tr_mat_mul_kernel");
}
else if (!d_mat.lhs().row_major() && result.row_major())
{
compressed_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_tr_mat_mul_kernel");
}
else
{
compressed_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle1()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_matrix_d_tr_mat_mul_kernel");
}
}
//
// triangular solves for compressed_matrix
//
template<typename NumericT>
__global__ void compressed_matrix_diagonal_kernel(
const unsigned int * row_indices,
const unsigned int * column_indices,
const NumericT * elements,
NumericT * result,
unsigned int size)
{
for (unsigned int row = blockDim.x * blockIdx.x + threadIdx.x;
row < size;
row += gridDim.x * blockDim.x)
{
NumericT diag = NumericT(0);
unsigned int row_end = row_indices[row+1];
for (unsigned int i = row_indices[row]; i < row_end; ++i)
{
unsigned int col_index = column_indices[i];
if (col_index == row)
{
diag = elements[i];
break;
}
}
result[row] = diag;
}
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const SparseMatrixT & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::unit_lower_tag)
{
csr_unit_lu_forward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_unit_lu_forward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const SparseMatrixT & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::lower_tag)
{
csr_lu_forward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_lu_forward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const SparseMatrixT & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::unit_upper_tag)
{
csr_unit_lu_backward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_unit_lu_backward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const SparseMatrixT & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::upper_tag)
{
csr_lu_backward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_lu_backward_kernel");
}
// transposed
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const matrix_expression<const SparseMatrixT, const SparseMatrixT, op_trans> & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::unit_lower_tag)
{
csr_trans_unit_lu_forward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.lhs().size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_trans_unit_lu_forward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const matrix_expression<const SparseMatrixT, const SparseMatrixT, op_trans> & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::lower_tag)
{
viennacl::vector<NumericT> diagonal(vec.size());
compressed_matrix_diagonal_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(diagonal),
static_cast<unsigned int>(mat.size1())
);
csr_trans_lu_forward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(diagonal),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.lhs().size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_trans_lu_forward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const matrix_expression<const SparseMatrixT, const SparseMatrixT, op_trans> & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::unit_upper_tag)
{
csr_trans_unit_lu_backward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.lhs().size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_trans_unit_lu_backward_kernel");
}
template<typename SparseMatrixT, typename NumericT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixT>::value>::type
inplace_solve(const matrix_expression<const SparseMatrixT, const SparseMatrixT, op_trans> & mat,
viennacl::vector_base<NumericT> & vec,
viennacl::linalg::upper_tag)
{
viennacl::vector<NumericT> diagonal(vec.size());
compressed_matrix_diagonal_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(diagonal),
static_cast<unsigned int>(mat.size1())
);
csr_trans_lu_backward_kernel<<<1, 128>>>(viennacl::cuda_arg<unsigned int>(mat.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(mat.lhs().handle2()),
viennacl::cuda_arg<NumericT>(mat.lhs().handle()),
viennacl::cuda_arg(diagonal),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(mat.lhs().size1())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("csr_trans_lu_backward_kernel");
}
namespace detail
{
//
// block solves
//
template<typename NumericT, unsigned int AlignmentV>
void block_inplace_solve(const matrix_expression<const compressed_matrix<NumericT, AlignmentV>,
const compressed_matrix<NumericT, AlignmentV>,
op_trans> & L,
viennacl::backend::mem_handle const & block_indices, vcl_size_t num_blocks,
vector_base<NumericT> const & /* L_diagonal */, //ignored
vector_base<NumericT> & vec,
viennacl::linalg::unit_lower_tag)
{
csr_block_trans_unit_lu_forward<<<num_blocks, 128>>>(viennacl::cuda_arg<unsigned int>(L.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(L.lhs().handle2()),
viennacl::cuda_arg<NumericT>(L.lhs().handle()),
viennacl::cuda_arg<unsigned int>(block_indices),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(L.lhs().size1())
);
}
template<typename NumericT, unsigned int AlignmentV>
void block_inplace_solve(const matrix_expression<const compressed_matrix<NumericT, AlignmentV>,
const compressed_matrix<NumericT, AlignmentV>,
op_trans> & U,
viennacl::backend::mem_handle const & block_indices, vcl_size_t num_blocks,
vector_base<NumericT> const & U_diagonal,
vector_base<NumericT> & vec,
viennacl::linalg::upper_tag)
{
csr_block_trans_lu_backward<<<num_blocks, 128>>>(viennacl::cuda_arg<unsigned int>(U.lhs().handle1()),
viennacl::cuda_arg<unsigned int>(U.lhs().handle2()),
viennacl::cuda_arg<NumericT>(U.lhs().handle()),
viennacl::cuda_arg(U_diagonal),
viennacl::cuda_arg<unsigned int>(block_indices),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(U.lhs().size1())
);
}
}
//
// Compressed Compressed Matrix
//
template<typename NumericT>
__global__ void compressed_compressed_matrix_vec_mul_kernel(
const unsigned int * row_jumper,
const unsigned int * row_indices,
const unsigned int * column_indices,
const NumericT * elements,
unsigned int nonzero_rows,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
unsigned int size_result,
NumericT beta)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < nonzero_rows;
i += gridDim.x * blockDim.x)
{
NumericT dot_prod = NumericT(0);
unsigned int row_end = row_jumper[i+1];
for (unsigned int j = row_jumper[i]; j < row_end; ++j)
dot_prod += elements[j] * x[column_indices[j] * inc_x + start_x];
unsigned int index = row_indices[i] * inc_result + start_result;
if (beta != 0) result[index] += alpha * dot_prod;
else result[index] = alpha * dot_prod;
}
}
template<typename NumericT>
void prod_impl(const viennacl::compressed_compressed_matrix<NumericT> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
if (beta < 0 || beta > 0)
viennacl::linalg::cuda::av(result, result, beta, 1, false, false);
else
result.clear();
compressed_compressed_matrix_vec_mul_kernel<<<128, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
static_cast<unsigned int>(mat.nnz1()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta
);
VIENNACL_CUDA_LAST_ERROR_CHECK("compressed_compressed_matrix_vec_mul_kernel");
}
//
// Coordinate Matrix
//
namespace detail
{
template<typename NumericT>
__global__ void coo_row_info_extractor( const unsigned int * coords, //(row_index, column_index)
const NumericT * elements,
const unsigned int * group_boundaries,
NumericT * result,
unsigned int option)
{
__shared__ unsigned int shared_rows[128];
__shared__ NumericT inter_results[128];
uint2 tmp;
NumericT val;
unsigned int last_index = blockDim.x - 1;
unsigned int group_start = group_boundaries[blockIdx.x];
unsigned int group_end = group_boundaries[blockIdx.x + 1];
unsigned int k_end = (group_end > group_start) ? 1 + (group_end - group_start - 1) / blockDim.x : 0; // -1 in order to have correct behavior if group_end - group_start == j * blockDim.x
unsigned int local_index = 0;
for (unsigned int k = 0; k < k_end; ++k)
{
local_index = group_start + k * blockDim.x + threadIdx.x;
tmp = (local_index < group_end) ? ((const uint2 *)coords)[local_index] : ::make_uint2(0, 0);
val = (local_index < group_end && (option != 3 || tmp.x == tmp.y) ) ? elements[local_index] : 0;
//check for carry from previous loop run:
if (threadIdx.x == 0 && k > 0)
{
if (tmp.x == shared_rows[last_index])
{
switch (option)
{
case 0: //inf-norm
case 3: //diagonal entry
val = max(val, fabs(inter_results[last_index]));
break;
case 1: //1-norm
val = fabs(val) + inter_results[last_index];
break;
case 2: //2-norm
val = sqrt(val * val + inter_results[last_index]);
break;
default:
break;
}
}
else
{
switch (option)
{
case 0: //inf-norm
case 1: //1-norm
case 3: //diagonal entry
result[shared_rows[last_index]] = inter_results[last_index];
break;
case 2: //2-norm
result[shared_rows[last_index]] = sqrt(inter_results[last_index]);
default:
break;
}
}
}
//segmented parallel reduction begin
__syncthreads();
shared_rows[threadIdx.x] = tmp.x;
switch (option)
{
case 0:
case 3:
inter_results[threadIdx.x] = val;
break;
case 1:
inter_results[threadIdx.x] = fabs(val);
break;
case 2:
inter_results[threadIdx.x] = val * val;
default:
break;
}
__syncthreads();
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
NumericT left = (threadIdx.x >= stride && tmp.x == shared_rows[threadIdx.x - stride]) ? inter_results[threadIdx.x - stride] : 0;
__syncthreads();
switch (option)
{
case 0: //inf-norm
case 3: //diagonal entry
inter_results[threadIdx.x] = max(inter_results[threadIdx.x], left);
break;
case 1: //1-norm
inter_results[threadIdx.x] += left;
break;
case 2: //2-norm
inter_results[threadIdx.x] += left;
break;
default:
break;
}
__syncthreads();
}
//segmented parallel reduction end
if (threadIdx.x != last_index &&
shared_rows[threadIdx.x] != shared_rows[threadIdx.x + 1] &&
inter_results[threadIdx.x] != 0)
{
result[tmp.x] = (option == 2) ? sqrt(inter_results[threadIdx.x]) : inter_results[threadIdx.x];
}
__syncthreads();
} //for k
if (local_index + 1 == group_end && inter_results[threadIdx.x] != 0)
result[tmp.x] = (option == 2) ? sqrt(inter_results[threadIdx.x]) : inter_results[threadIdx.x];
}
template<typename NumericT, unsigned int AlignmentV>
void row_info(coordinate_matrix<NumericT, AlignmentV> const & mat,
vector_base<NumericT> & vec,
viennacl::linalg::detail::row_info_types info_selector)
{
coo_row_info_extractor<<<64, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle12()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(info_selector)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coo_row_info_extractor");
}
} //namespace detail
template<typename NumericT>
__global__ void coordinate_matrix_vec_mul_kernel(const unsigned int * coords, //(row_index, column_index)
const NumericT * elements,
const unsigned int * group_boundaries,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
NumericT beta)
{
__shared__ unsigned int shared_rows[128];
__shared__ NumericT inter_results[128];
uint2 tmp;
NumericT val;
unsigned int group_start = group_boundaries[blockIdx.x];
unsigned int group_end = group_boundaries[blockIdx.x + 1];
unsigned int k_end = (group_end > group_start) ? 1 + (group_end - group_start - 1) / blockDim.x : 0; // -1 in order to have correct behavior if group_end - group_start == j * blockDim.x
unsigned int local_index = 0;
for (unsigned int k = 0; k < k_end; ++k)
{
local_index = group_start + k * blockDim.x + threadIdx.x;
tmp = (local_index < group_end) ? ((const uint2 *)coords)[local_index] : ::make_uint2(0, 0);
val = (local_index < group_end) ? elements[local_index] * x[tmp.y * inc_x + start_x] : 0;
//check for carry from previous loop run:
if (threadIdx.x == 0 && k > 0)
{
if (tmp.x == shared_rows[blockDim.x-1])
val += inter_results[blockDim.x-1];
else if (beta != 0)
result[shared_rows[blockDim.x-1] * inc_result + start_result] += alpha * inter_results[blockDim.x-1];
else
result[shared_rows[blockDim.x-1] * inc_result + start_result] = alpha * inter_results[blockDim.x-1];
}
//segmented parallel reduction begin
__syncthreads();
shared_rows[threadIdx.x] = tmp.x;
inter_results[threadIdx.x] = val;
NumericT left = 0;
__syncthreads();
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
left = (threadIdx.x >= stride && tmp.x == shared_rows[threadIdx.x - stride]) ? inter_results[threadIdx.x - stride] : 0;
__syncthreads();
inter_results[threadIdx.x] += left;
__syncthreads();
}
//segmented parallel reduction end
if (local_index < group_end - 1 && threadIdx.x < blockDim.x-1 &&
shared_rows[threadIdx.x] != shared_rows[threadIdx.x + 1])
{
if (beta != 0) result[tmp.x * inc_result + start_result] += alpha * inter_results[threadIdx.x];
else result[tmp.x * inc_result + start_result] = alpha * inter_results[threadIdx.x];
}
__syncthreads();
} //for k
if (local_index + 1 == group_end) {
if (beta != 0) result[tmp.x * inc_result + start_result] += alpha * inter_results[threadIdx.x];
else result[tmp.x * inc_result + start_result] = alpha * inter_results[threadIdx.x];
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::coordinate_matrix<NumericT, AlignmentV> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
if (beta < 0 || beta > 0)
viennacl::linalg::cuda::av(result, result, beta, 1, false, false);
else
result.clear();
coordinate_matrix_vec_mul_kernel<<<64, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle12()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
beta
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_vec_mul_kernel");
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void coordinate_matrix_d_mat_mul_kernel(const unsigned int * coords, //(row_index, column_index)
const NumericT * elements,
const unsigned int * group_boundaries,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
__shared__ unsigned int shared_rows[128];
__shared__ NumericT inter_results[128];
uint2 tmp;
NumericT val;
unsigned int group_start = group_boundaries[blockIdx.x];
unsigned int group_end = group_boundaries[blockIdx.x + 1];
unsigned int k_end = (group_end > group_start) ? 1 + (group_end - group_start - 1) / blockDim.x : 0; // -1 in order to have correct behavior if group_end - group_start == j * blockDim.x
unsigned int local_index = 0;
for (unsigned int result_col = 0; result_col < result_col_size; ++result_col)
{
for (unsigned int k = 0; k < k_end; ++k)
{
local_index = group_start + k * blockDim.x + threadIdx.x;
tmp = (local_index < group_end) ? ((const uint2 *)coords)[local_index] : ::make_uint2(0, 0);
val = (local_index < group_end) ? elements[local_index] * d_mat[DMatIndexT::apply(tmp.y, result_col,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols) ] : 0;
//check for carry from previous loop run:
if (threadIdx.x == 0 && k > 0)
{
if (tmp.x == shared_rows[blockDim.x-1])
val += inter_results[blockDim.x-1];
else
result[ResultIndexT::apply(shared_rows[blockDim.x-1], result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = inter_results[blockDim.x-1];
}
//segmented parallel reduction begin
__syncthreads();
shared_rows[threadIdx.x] = tmp.x;
inter_results[threadIdx.x] = val;
NumericT left = 0;
__syncthreads();
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
left = (threadIdx.x >= stride && tmp.x == shared_rows[threadIdx.x - stride]) ? inter_results[threadIdx.x - stride] : 0;
__syncthreads();
inter_results[threadIdx.x] += left;
__syncthreads();
}
//segmented parallel reduction end
if (local_index < group_end && threadIdx.x < blockDim.x-1 &&
shared_rows[threadIdx.x] != shared_rows[threadIdx.x + 1])
{
result[ResultIndexT::apply(tmp.x, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = inter_results[threadIdx.x];
}
__syncthreads();
} //for k
if (local_index + 1 == group_end)
result[ResultIndexT::apply(tmp.x, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = inter_results[threadIdx.x];
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::coordinate_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_base<NumericT> & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.row_major() && result.row_major())
{
coordinate_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_mat_mul_kernel");
}
else if (d_mat.row_major() && !result.row_major())
{
coordinate_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_mat_mul_kernel");
}
else if (!d_mat.row_major() && result.row_major())
{
coordinate_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_mat_mul_kernel");
}
else
{
coordinate_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_mat_mul_kernel");
}
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void coordinate_matrix_d_tr_mat_mul_kernel(const unsigned int * coords, //(row_index, column_index)
const NumericT * elements,
const unsigned int * group_boundaries,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
__shared__ unsigned int shared_rows[128];
__shared__ NumericT inter_results[128];
uint2 tmp;
NumericT val;
unsigned int group_start = group_boundaries[blockIdx.x];
unsigned int group_end = group_boundaries[blockIdx.x + 1];
unsigned int k_end = (group_end > group_start) ? 1 + (group_end - group_start - 1) / blockDim.x : 0; // -1 in order to have correct behavior if group_end - group_start == j * blockDim.x
unsigned int local_index = 0;
for (unsigned int result_col = 0; result_col < result_col_size; ++result_col)
{
for (unsigned int k = 0; k < k_end; ++k)
{
local_index = group_start + k * blockDim.x + threadIdx.x;
tmp = (local_index < group_end) ? ((const uint2 *)coords)[local_index] : ::make_uint2(0, 0);
val = (local_index < group_end) ? elements[local_index] * d_mat[DMatIndexT::apply(result_col, tmp.y,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols)] : 0;
//check for carry from previous loop run:
if (threadIdx.x == 0 && k > 0)
{
if (tmp.x == shared_rows[blockDim.x-1])
val += inter_results[blockDim.x-1];
else
result[ResultIndexT::apply(shared_rows[blockDim.x-1], result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = inter_results[blockDim.x-1];
}
//segmented parallel reduction begin
__syncthreads();
shared_rows[threadIdx.x] = tmp.x;
inter_results[threadIdx.x] = val;
NumericT left = 0;
__syncthreads();
for (unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
left = (threadIdx.x >= stride && tmp.x == shared_rows[threadIdx.x - stride]) ? inter_results[threadIdx.x - stride] : 0;
__syncthreads();
inter_results[threadIdx.x] += left;
__syncthreads();
}
//segmented parallel reduction end
if (local_index < group_end && threadIdx.x < blockDim.x-1 &&
shared_rows[threadIdx.x] != shared_rows[threadIdx.x + 1])
{
result[ ResultIndexT::apply(tmp.x, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = inter_results[threadIdx.x];
}
__syncthreads();
} //for k
if (local_index + 1 == group_end)
result[ ResultIndexT::apply(tmp.x, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = inter_results[threadIdx.x];
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::coordinate_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
const viennacl::matrix_base<NumericT>,
viennacl::op_trans > & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.lhs().row_major() && result.row_major())
{
coordinate_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_tr_mat_mul_kernel");
}
else if (d_mat.lhs().row_major() && !result.row_major())
{
coordinate_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_tr_mat_mul_kernel");
}
else if (!d_mat.lhs().row_major() && result.row_major())
{
coordinate_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_tr_mat_mul_kernel");
}
else
{
coordinate_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<64, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle12()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
viennacl::cuda_arg<unsigned int>(sp_mat.handle3()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("coordinate_matrix_d_tr_mat_mul_kernel");
}
}
//
// ELL Matrix
//
template<typename AlphaBetaHandlerT, typename NumericT>
__global__ void ell_matrix_vec_mul_kernel(const unsigned int * coords,
const NumericT * elements,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
NumericT beta,
unsigned int row_num,
unsigned int col_num,
unsigned int internal_row_num,
unsigned int items_per_row,
unsigned int aligned_items_per_row
)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int row_id = glb_id; row_id < row_num; row_id += glb_sz)
{
NumericT sum = 0;
unsigned int offset = row_id;
for (unsigned int item_id = 0; item_id < items_per_row; item_id++, offset += internal_row_num)
{
NumericT val = elements[offset];
if (val != NumericT(0))
{
int col = coords[offset];
sum += x[col * inc_x + start_x] * val;
}
}
AlphaBetaHandlerT::apply(result[row_id * inc_result + start_result], alpha, sum, beta);
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::ell_matrix<NumericT, AlignmentV> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
ell_matrix_vec_mul_kernel<detail::spmv_alpha_beta><<<256, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
beta,
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.size2()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.maxnnz()),
static_cast<unsigned int>(mat.internal_maxnnz())
);
else
ell_matrix_vec_mul_kernel<detail::spmv_pure><<<256, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
beta,
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.size2()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.maxnnz()),
static_cast<unsigned int>(mat.internal_maxnnz())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_vec_mul_kernel");
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void ell_matrix_d_mat_mul_kernel(const unsigned int * sp_mat_coords,
const NumericT * sp_mat_elements,
unsigned int sp_mat_row_num,
unsigned int sp_mat_col_num,
unsigned int sp_mat_internal_row_num,
unsigned int sp_mat_items_per_row,
unsigned int sp_mat_aligned_items_per_row,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for ( unsigned int rc = glb_id; rc < (sp_mat_row_num * d_mat_col_size); rc += glb_sz)
{
unsigned int row = rc % sp_mat_row_num;
unsigned int col = rc / sp_mat_row_num;
unsigned int offset = row;
NumericT r = (NumericT)0;
for (unsigned int k = 0; k < sp_mat_items_per_row; k++, offset += sp_mat_internal_row_num)
{
unsigned int j = sp_mat_coords[offset];
NumericT x = static_cast<NumericT>(sp_mat_elements[offset]);
if (x != (NumericT)0)
{
NumericT y = d_mat[ DMatIndexT::apply(j, col,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols) ];
r += x*y;
}
}
result [ ResultIndexT::apply(row, col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = r;
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::ell_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_base<NumericT> & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.row_major() && result.row_major())
{
ell_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_mat_mul_kernel");
}
else if (d_mat.row_major() && !result.row_major())
{
ell_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_mat_mul_kernel");
}
else if (!d_mat.row_major() && result.row_major())
{
ell_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_mat_mul_kernel");
}
else
{
ell_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_mat_mul_kernel");
}
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT >
__global__ void ell_matrix_d_tr_mat_mul_kernel(const unsigned int * sp_mat_coords,
const NumericT * sp_mat_elements,
unsigned int sp_mat_row_num,
unsigned int sp_mat_col_num,
unsigned int sp_mat_internal_row_num,
unsigned int sp_mat_items_per_row,
unsigned int sp_mat_aligned_items_per_row,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for ( unsigned int rc = glb_id; rc < (sp_mat_row_num * d_mat_row_size); rc += glb_sz)
{
unsigned int row = rc % sp_mat_row_num;
unsigned int col = rc / sp_mat_row_num;
unsigned int offset = row;
NumericT r = (NumericT)0;
for (unsigned int k = 0; k < sp_mat_items_per_row; k++, offset += sp_mat_internal_row_num)
{
unsigned int j = sp_mat_coords[offset];
NumericT x = static_cast<NumericT>(sp_mat_elements[offset]);
if (x != (NumericT)0)
{
NumericT y = d_mat[ DMatIndexT::apply(col, j,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols) ];
r += x*y;
}
}
result [ ResultIndexT::apply(row, col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols) ] = r;
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::ell_matrix<NumericT, AlignmentV> & sp_mat,
const viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
const viennacl::matrix_base<NumericT>,
viennacl::op_trans > & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.lhs().row_major() && result.row_major())
{
ell_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_tr_mat_mul_kernel");
}
else if (d_mat.lhs().row_major() && !result.row_major())
{
ell_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_tr_mat_mul_kernel");
}
else if (!d_mat.lhs().row_major() && result.row_major())
{
ell_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_tr_mat_mul_kernel");
}
else
{
ell_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<128, 128>>>
(viennacl::cuda_arg<unsigned int>(sp_mat.handle2()),
viennacl::cuda_arg<NumericT>(sp_mat.handle()),
static_cast<unsigned int>(sp_mat.size1()),
static_cast<unsigned int>(sp_mat.size2()),
static_cast<unsigned int>(sp_mat.internal_size1()),
static_cast<unsigned int>(sp_mat.maxnnz()),
static_cast<unsigned int>(sp_mat.internal_maxnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("ell_matrix_d_tr_mat_mul_kernel");
}
}
//
// SELL-C-\sigma Matrix
//
template<typename AlphaBetaHandlerT, typename NumericT>
__global__ void sliced_ell_matrix_vec_mul_kernel(const unsigned int * columns_per_block,
const unsigned int * column_indices,
const unsigned int * block_start,
const NumericT * elements,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
unsigned int size_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
unsigned int size_result,
NumericT beta,
unsigned int block_size)
{
unsigned int blocks_per_threadblock = blockDim.x / block_size;
unsigned int id_in_block = threadIdx.x % block_size;
unsigned int num_blocks = (size_result - 1) / block_size + 1;
unsigned int global_warp_count = blocks_per_threadblock * gridDim.x;
unsigned int global_warp_id = blocks_per_threadblock * blockIdx.x + threadIdx.x / block_size;
for (unsigned int block_idx = global_warp_id; block_idx < num_blocks; block_idx += global_warp_count)
{
unsigned int row = block_idx * block_size + id_in_block;
unsigned int offset = block_start[block_idx];
unsigned int num_columns = columns_per_block[block_idx];
NumericT sum = 0;
for (unsigned int item_id = 0; item_id < num_columns; item_id++)
{
unsigned int index = offset + item_id * block_size + id_in_block;
NumericT val = elements[index];
sum += val ? (x[column_indices[index] * inc_x + start_x] * val) : 0;
}
if (row < size_result)
AlphaBetaHandlerT::apply(result[row * inc_result + start_result], alpha, sum, beta);
}
}
template<typename NumericT, typename IndexT>
void prod_impl(const viennacl::sliced_ell_matrix<NumericT, IndexT> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
sliced_ell_matrix_vec_mul_kernel<detail::spmv_alpha_beta><<<256, 256>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
static_cast<unsigned int>(vec.size()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta,
static_cast<unsigned int>(mat.rows_per_block())
);
else
sliced_ell_matrix_vec_mul_kernel<detail::spmv_pure><<<256, 256>>>(viennacl::cuda_arg<unsigned int>(mat.handle1()),
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
static_cast<unsigned int>(vec.size()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
static_cast<unsigned int>(result.size()),
beta,
static_cast<unsigned int>(mat.rows_per_block())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("sliced_ell_matrix_vec_mul_kernel");
}
//
// Hybrid Matrix
//
template<typename AlphaBetaHandlerT, typename NumericT>
__global__ void hyb_matrix_vec_mul_kernel(const unsigned int * ell_coords,
const NumericT * ell_elements,
const unsigned int * csr_rows,
const unsigned int * csr_cols,
const NumericT * csr_elements,
const NumericT * x,
unsigned int start_x,
unsigned int inc_x,
NumericT alpha,
NumericT * result,
unsigned int start_result,
unsigned int inc_result,
NumericT beta,
unsigned int row_num,
unsigned int internal_row_num,
unsigned int items_per_row,
unsigned int aligned_items_per_row
)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int row_id = glb_id; row_id < row_num; row_id += glb_sz)
{
NumericT sum = 0;
unsigned int offset = row_id;
for (unsigned int item_id = 0; item_id < items_per_row; item_id++, offset += internal_row_num)
{
NumericT val = ell_elements[offset];
if (val != NumericT(0))
{
int col = ell_coords[offset];
sum += (x[col * inc_x + start_x] * val);
}
}
unsigned int col_begin = csr_rows[row_id];
unsigned int col_end = csr_rows[row_id + 1];
for (unsigned int item_id = col_begin; item_id < col_end; item_id++)
sum += x[csr_cols[item_id] * inc_x + start_x] * csr_elements[item_id];
AlphaBetaHandlerT::apply(result[row_id * inc_result + start_result], alpha, sum, beta);
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::hyb_matrix<NumericT, AlignmentV> & mat,
const viennacl::vector_base<NumericT> & vec,
NumericT alpha,
viennacl::vector_base<NumericT> & result,
NumericT beta)
{
if (alpha < NumericT(1) || alpha > NumericT(1) || beta < 0 || beta > 0)
hyb_matrix_vec_mul_kernel<detail::spmv_alpha_beta><<<256, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
beta,
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz())
);
else
hyb_matrix_vec_mul_kernel<detail::spmv_pure><<<256, 128>>>(viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(vec.start()),
static_cast<unsigned int>(vec.stride()),
alpha,
viennacl::cuda_arg(result),
static_cast<unsigned int>(result.start()),
static_cast<unsigned int>(result.stride()),
beta,
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void hyb_matrix_d_mat_mul_kernel(const unsigned int * ell_coords,
const NumericT * ell_elements,
const unsigned int * csr_rows,
const unsigned int * csr_cols,
const NumericT * csr_elements,
unsigned int row_num,
unsigned int internal_row_num,
unsigned int items_per_row,
unsigned int aligned_items_per_row,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int result_col = 0; result_col < result_col_size; ++result_col)
{
for (unsigned int row_id = glb_id; row_id < row_num; row_id += glb_sz)
{
NumericT sum = 0;
unsigned int offset = row_id;
for (unsigned int item_id = 0; item_id < items_per_row; item_id++, offset += internal_row_num)
{
NumericT val = ell_elements[offset];
if (val != 0.0f)
{
sum += d_mat[DMatIndexT::apply(ell_coords[offset], result_col,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols)] * val;
}
}
unsigned int col_begin = csr_rows[row_id];
unsigned int col_end = csr_rows[row_id + 1];
for (unsigned int item_id = col_begin; item_id < col_end; item_id++)
{
sum += d_mat[DMatIndexT::apply(csr_cols[item_id], result_col,
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols)] * csr_elements[item_id];
}
result[ResultIndexT::apply(row_id, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = sum;
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::hyb_matrix<NumericT, AlignmentV> & mat,
const viennacl::matrix_base<NumericT> & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.row_major() && result.row_major())
{
hyb_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else if (d_mat.row_major() && !result.row_major())
{
hyb_matrix_d_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else if (!d_mat.row_major() && result.row_major())
{
hyb_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else
{
hyb_matrix_d_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat),
static_cast<unsigned int>(viennacl::traits::start1(d_mat)), static_cast<unsigned int>(viennacl::traits::start2(d_mat)),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat)), static_cast<unsigned int>(viennacl::traits::stride2(d_mat)),
static_cast<unsigned int>(viennacl::traits::size1(d_mat)), static_cast<unsigned int>(viennacl::traits::size2(d_mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat)),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
}
template<typename DMatIndexT, typename ResultIndexT, typename NumericT>
__global__ void hyb_matrix_d_tr_mat_mul_kernel(const unsigned int * ell_coords,
const NumericT * ell_elements,
const unsigned int * csr_rows,
const unsigned int * csr_cols,
const NumericT * csr_elements,
unsigned int row_num,
unsigned int internal_row_num,
unsigned int items_per_row,
unsigned int aligned_items_per_row,
const NumericT * d_mat,
unsigned int d_mat_row_start,
unsigned int d_mat_col_start,
unsigned int d_mat_row_inc,
unsigned int d_mat_col_inc,
unsigned int d_mat_row_size,
unsigned int d_mat_col_size,
unsigned int d_mat_internal_rows,
unsigned int d_mat_internal_cols,
NumericT * result,
unsigned int result_row_start,
unsigned int result_col_start,
unsigned int result_row_inc,
unsigned int result_col_inc,
unsigned int result_row_size,
unsigned int result_col_size,
unsigned int result_internal_rows,
unsigned int result_internal_cols)
{
unsigned int glb_id = blockDim.x * blockIdx.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int result_col = 0; result_col < result_col_size; ++result_col)
{
for (unsigned int row_id = glb_id; row_id < row_num; row_id += glb_sz)
{
NumericT sum = 0;
unsigned int offset = row_id;
for (unsigned int item_id = 0; item_id < items_per_row; item_id++, offset += internal_row_num)
{
NumericT val = ell_elements[offset];
if (val != 0.0f)
{
sum += d_mat[DMatIndexT::apply(result_col, ell_coords[offset],
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols)] * val;
}
}
unsigned int col_begin = csr_rows[row_id];
unsigned int col_end = csr_rows[row_id + 1];
for (unsigned int item_id = col_begin; item_id < col_end; item_id++)
{
sum += d_mat[DMatIndexT::apply(result_col, csr_cols[item_id],
d_mat_row_start, d_mat_row_inc,
d_mat_col_start, d_mat_col_inc,
d_mat_internal_rows, d_mat_internal_cols)] * csr_elements[item_id];
}
result[ResultIndexT::apply(row_id, result_col,
result_row_start, result_row_inc,
result_col_start, result_col_inc,
result_internal_rows, result_internal_cols)] = sum;
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void prod_impl(const viennacl::hyb_matrix<NumericT, AlignmentV> & mat,
const viennacl::matrix_expression< const viennacl::matrix_base<NumericT>,
const viennacl::matrix_base<NumericT>,
viennacl::op_trans > & d_mat,
viennacl::matrix_base<NumericT> & result)
{
if (d_mat.lhs().row_major() && result.row_major())
{
hyb_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<row_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else if (d_mat.lhs().row_major() && !result.row_major())
{
hyb_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<row_major>, mat_mult_matrix_index<column_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else if (!d_mat.lhs().row_major() && result.row_major())
{
hyb_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<row_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
else
{
hyb_matrix_d_tr_mat_mul_kernel<mat_mult_matrix_index<column_major>, mat_mult_matrix_index<column_major> ><<<256, 128>>>(
viennacl::cuda_arg<unsigned int>(mat.handle2()),
viennacl::cuda_arg<NumericT>(mat.handle()),
viennacl::cuda_arg<unsigned int>(mat.handle3()),
viennacl::cuda_arg<unsigned int>(mat.handle4()),
viennacl::cuda_arg<NumericT>(mat.handle5()),
static_cast<unsigned int>(mat.size1()),
static_cast<unsigned int>(mat.internal_size1()),
static_cast<unsigned int>(mat.ell_nnz()),
static_cast<unsigned int>(mat.internal_ellnnz()),
viennacl::cuda_arg(d_mat.lhs()),
static_cast<unsigned int>(viennacl::traits::start1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::start2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::stride1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::stride2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::size2(d_mat.lhs())),
static_cast<unsigned int>(viennacl::traits::internal_size1(d_mat.lhs())), static_cast<unsigned int>(viennacl::traits::internal_size2(d_mat.lhs())),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start1(result)), static_cast<unsigned int>(viennacl::traits::start2(result)),
static_cast<unsigned int>(viennacl::traits::stride1(result)), static_cast<unsigned int>(viennacl::traits::stride2(result)),
static_cast<unsigned int>(viennacl::traits::size1(result)), static_cast<unsigned int>(viennacl::traits::size2(result)),
static_cast<unsigned int>(viennacl::traits::internal_size1(result)), static_cast<unsigned int>(viennacl::traits::internal_size2(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("hyb_matrix_vec_mul_kernel");
}
}
} // namespace cuda
} //namespace linalg
} //namespace viennacl
#endif
viennacl::hyb_matrix::internal_ellnnz
vcl_size_t internal_ellnnz() const
Definition: hyb_matrix.hpp:101
viennacl::hyb_matrix
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
Definition: forwards.h:406
viennacl::enable_if
Simple enable-if variant that uses the SFINAE pattern.
Definition: enable_if.hpp:30
viennacl::sliced_ell_matrix::handle2
handle_type & handle2()
Definition: sliced_ell_matrix.hpp:112
viennacl::ell_matrix::size1
vcl_size_t size1() const
Definition: ell_matrix.hpp:91
viennacl::linalg::cuda::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:253
viennacl::ell_matrix::handle2
handle_type & handle2()
Definition: ell_matrix.hpp:103
viennacl::traits::stride1
result_of::size_type< matrix_base< NumericT > >::type stride1(matrix_base< NumericT > const &s)
Definition: stride.hpp:55
viennacl::linalg::sum
viennacl::scalar_expression< const viennacl::vector_base< NumericT >, const viennacl::vector_base< NumericT >, viennacl::op_sum > sum(viennacl::vector_base< NumericT > const &x)
User interface function for computing the sum of all elements of a vector.
Definition: sum.hpp:45
viennacl::linalg::detail::row_info_types
row_info_types
Definition: forwards.h:837
viennacl::linalg::cuda::ell_matrix_vec_mul_kernel
__global__ void ell_matrix_vec_mul_kernel(const unsigned int *coords, const NumericT *elements, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, NumericT beta, unsigned int row_num, unsigned int col_num, unsigned int internal_row_num, unsigned int items_per_row, unsigned int aligned_items_per_row)
Definition: sparse_matrix_operations.hpp:1747
viennacl::hyb_matrix::handle3
const handle_type & handle3() const
Definition: hyb_matrix.hpp:107
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
viennacl::compressed_compressed_matrix::handle2
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_compressed_matrix.hpp:471
viennacl::compressed_compressed_matrix::handle1
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_compressed_matrix.hpp:469
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::hyb_matrix::handle
const handle_type & handle() const
Definition: hyb_matrix.hpp:105
viennacl::traits::internal_size1
vcl_size_t internal_size1(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per row of a ViennaCL matrix...
Definition: size.hpp:386
viennacl::coordinate_matrix::handle12
const handle_type & handle12() const
Returns the OpenCL handle to the (row, column) index array.
Definition: coordinate_matrix.hpp:366
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
viennacl::linalg::cuda::coordinate_matrix_d_mat_mul_kernel
__global__ void coordinate_matrix_d_mat_mul_kernel(const unsigned int *coords, const NumericT *elements, const unsigned int *group_boundaries, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:1352
viennacl::linalg::cuda::av
void av(vector_base< NumericT > &vec1, vector_base< NumericT > const &vec2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: vector_operations.hpp:144
viennacl::linalg::cuda::detail::spmv_pure::apply
static __device__ void apply(NumericT &result, NumericT alpha, NumericT Ax, NumericT beta)
Definition: sparse_matrix_operations.hpp:124
viennacl::hyb_matrix::internal_size1
vcl_size_t internal_size1() const
Definition: hyb_matrix.hpp:95
viennacl::traits::internal_size2
vcl_size_t internal_size2(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per column of a ViennaCL matrix...
Definition: size.hpp:394
viennacl::linalg::cuda::hyb_matrix_vec_mul_kernel
__global__ void hyb_matrix_vec_mul_kernel(const unsigned int *ell_coords, const NumericT *ell_elements, const unsigned int *csr_rows, const unsigned int *csr_cols, const NumericT *csr_elements, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, NumericT beta, unsigned int row_num, unsigned int internal_row_num, unsigned int items_per_row, unsigned int aligned_items_per_row)
Definition: sparse_matrix_operations.hpp:2298
viennacl::matrix_base< NumericT >
viennacl::matrix_expression
Expression template class for representing a tree of expressions which ultimately result in a matrix...
Definition: forwards.h:341
viennacl::traits::stride
result_of::size_type< viennacl::vector_base< T > >::type stride(viennacl::vector_base< T > const &s)
Definition: stride.hpp:45
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::ell_matrix::size2
vcl_size_t size2() const
Definition: ell_matrix.hpp:92
viennacl::traits::start1
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:65
viennacl::linalg::cuda::detail::row_info
void row_info(compressed_matrix< NumericT, AligmentV > const &mat, vector_base< NumericT > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:107
viennacl::linalg::detail::spai::dot_prod
void dot_prod(MatrixT const &A, unsigned int beg_ind, NumericT &res)
Dot prod of particular column of martix A with it's self starting at a certain index beg_ind...
Definition: qr.hpp:182
viennacl::hyb_matrix::handle4
const handle_type & handle4() const
Definition: hyb_matrix.hpp:108
viennacl::linalg::cuda::coordinate_matrix_vec_mul_kernel
__global__ void coordinate_matrix_vec_mul_kernel(const unsigned int *coords, const NumericT *elements, const unsigned int *group_boundaries, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, NumericT beta)
Definition: sparse_matrix_operations.hpp:1239
viennacl::sliced_ell_matrix::rows_per_block
vcl_size_t rows_per_block() const
Definition: sliced_ell_matrix.hpp:104
viennacl::linalg::cuda::prod_impl
void prod_impl(const matrix_base< NumericT > &mat, bool mat_transpose, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:1464
viennacl::linalg::cuda::compressed_matrix_diagonal_kernel
__global__ void compressed_matrix_diagonal_kernel(const unsigned int *row_indices, const unsigned int *column_indices, const NumericT *elements, NumericT *result, unsigned int size)
Definition: sparse_matrix_operations.hpp:746
viennacl::linalg::cuda::detail::spmv_alpha_beta::apply
static __device__ void apply(NumericT &result, NumericT alpha, NumericT Ax, NumericT beta)
Definition: sparse_matrix_operations.hpp:130
viennacl::traits::size2
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:201
viennacl::matrix_expression::size1
vcl_size_t size1() const
Returns the size of the result vector.
Definition: matrix.hpp:72
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
viennacl::linalg::cuda::mat_mult_matrix_index
Helper struct for accessing an element of a row- or column-major matrix.
Definition: sparse_matrix_operations.hpp:403
viennacl::ell_matrix::internal_size1
vcl_size_t internal_size1() const
Definition: ell_matrix.hpp:88
viennacl::compressed_matrix::nnz
const vcl_size_t & nnz() const
Returns the number of nonzero entries.
Definition: compressed_matrix.hpp:931
viennacl::sliced_ell_matrix::handle
handle_type & handle()
Definition: sliced_ell_matrix.hpp:118
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::cuda::detail::csr_row_info_extractor_kernel
__global__ void csr_row_info_extractor_kernel(const unsigned int *row_indices, const unsigned int *column_indices, const NumericT *elements, NumericT *result, unsigned int size, unsigned int option)
Definition: sparse_matrix_operations.hpp:54
viennacl::vector_base::stride
size_type stride() const
Returns the stride within the buffer (in multiples of sizeof(NumericT))
Definition: vector_def.hpp:124
viennacl::hyb_matrix::handle2
const handle_type & handle2() const
Definition: hyb_matrix.hpp:106
viennacl::compressed_compressed_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_compressed_matrix.hpp:475
viennacl::linalg::cuda::detail::coo_row_info_extractor
__global__ void coo_row_info_extractor(const unsigned int *coords, const NumericT *elements, const unsigned int *group_boundaries, NumericT *result, unsigned int option)
Definition: sparse_matrix_operations.hpp:1095
viennacl::linalg::cuda::hyb_matrix_d_tr_mat_mul_kernel
__global__ void hyb_matrix_d_tr_mat_mul_kernel(const unsigned int *ell_coords, const NumericT *ell_elements, const unsigned int *csr_rows, const unsigned int *csr_cols, const NumericT *csr_elements, unsigned int row_num, unsigned int internal_row_num, unsigned int items_per_row, unsigned int aligned_items_per_row, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:2606
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::traits::start2
result_of::size_type< T >::type start2(T const &obj)
Definition: start.hpp:84
viennacl::ell_matrix
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
viennacl::linalg::cuda::compressed_compressed_matrix_vec_mul_kernel
__global__ void compressed_compressed_matrix_vec_mul_kernel(const unsigned int *row_jumper, const unsigned int *row_indices, const unsigned int *column_indices, const NumericT *elements, unsigned int nonzero_rows, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, unsigned int size_result, NumericT beta)
Definition: sparse_matrix_operations.hpp:1016
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::sliced_ell_matrix
Sparse matrix class using the sliced ELLPACK with parameters C, .
Definition: forwards.h:403
viennacl::linalg::cuda::ell_matrix_d_mat_mul_kernel
__global__ void ell_matrix_d_mat_mul_kernel(const unsigned int *sp_mat_coords, const NumericT *sp_mat_elements, unsigned int sp_mat_row_num, unsigned int sp_mat_col_num, unsigned int sp_mat_internal_row_num, unsigned int sp_mat_items_per_row, unsigned int sp_mat_aligned_items_per_row, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:1841
viennacl::compressed_compressed_matrix::handle3
const handle_type & handle3() const
Returns the OpenCL handle to the row index array.
Definition: compressed_compressed_matrix.hpp:473
viennacl::compressed_compressed_matrix
A sparse square matrix in compressed sparse rows format optimized for the case that only a few rows c...
Definition: compressed_compressed_matrix.hpp:265
viennacl::compressed_matrix::handle2
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:938
viennacl::linalg::cuda::hyb_matrix_d_mat_mul_kernel
__global__ void hyb_matrix_d_mat_mul_kernel(const unsigned int *ell_coords, const NumericT *ell_elements, const unsigned int *csr_rows, const unsigned int *csr_cols, const NumericT *csr_elements, unsigned int row_num, unsigned int internal_row_num, unsigned int items_per_row, unsigned int aligned_items_per_row, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:2408
viennacl::vector_base< NumericT >
viennacl::linalg::cuda::compressed_matrix_vec_mul_kernel
__global__ void compressed_matrix_vec_mul_kernel(const unsigned int *row_indices, const unsigned int *column_indices, const NumericT *elements, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, unsigned int size_result, NumericT beta)
Definition: sparse_matrix_operations.hpp:138
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::vector< NumericT >
viennacl::diag
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
Definition: matrix.hpp:895
viennacl::linalg::cuda::mat_mult_matrix_index::apply
static __device__ unsigned int apply(unsigned int i, unsigned int j, unsigned int row_start, unsigned int row_inc, unsigned int col_start, unsigned int col_inc, unsigned int internal_rows, unsigned int internal_cols)
Definition: sparse_matrix_operations.hpp:405
viennacl::ell_matrix::maxnnz
vcl_size_t maxnnz() const
Definition: ell_matrix.hpp:95
viennacl::linalg::cuda::sliced_ell_matrix_vec_mul_kernel
__global__ void sliced_ell_matrix_vec_mul_kernel(const unsigned int *columns_per_block, const unsigned int *column_indices, const unsigned int *block_start, const NumericT *elements, const NumericT *x, unsigned int start_x, unsigned int inc_x, unsigned int size_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, unsigned int size_result, NumericT beta, unsigned int block_size)
Definition: sparse_matrix_operations.hpp:2197
viennacl::linalg::cuda::detail::block_inplace_solve
void block_inplace_solve(const matrix_expression< const compressed_matrix< NumericT, AlignmentV >, const compressed_matrix< NumericT, AlignmentV >, op_trans > &L, viennacl::backend::mem_handle const &block_indices, vcl_size_t num_blocks, vector_base< NumericT > const &, vector_base< NumericT > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:970
viennacl::traits::stride2
result_of::size_type< matrix_base< NumericT > >::type stride2(matrix_base< NumericT > const &s)
Definition: stride.hpp:65
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::coordinate_matrix::handle3
const handle_type & handle3() const
Returns the OpenCL handle to the group start index array.
Definition: coordinate_matrix.hpp:370
viennacl::linalg::cuda::detail::spmv_pure
Definition: sparse_matrix_operations.hpp:121
sparse_matrix_operations_solve.hpp
Implementations of direct triangular solvers for sparse matrices using CUDA.
viennacl::ell_matrix::handle
handle_type & handle()
Definition: ell_matrix.hpp:100
viennacl::compressed_matrix::handle3
const handle_type & handle3() const
Returns the OpenCL handle to the row block array.
Definition: compressed_matrix.hpp:940
spgemm_rmerge.hpp
viennacl::vector_base::clear
void clear()
Resets all entries to zero. Does not change the size of the vector.
Definition: vector.hpp:875
common.hpp
Common routines for CUDA execution.
viennacl::coordinate_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: coordinate_matrix.hpp:368
viennacl::linalg::cuda::ell_matrix_d_tr_mat_mul_kernel
__global__ void ell_matrix_d_tr_mat_mul_kernel(const unsigned int *sp_mat_coords, const NumericT *sp_mat_elements, unsigned int sp_mat_row_num, unsigned int sp_mat_col_num, unsigned int sp_mat_internal_row_num, unsigned int sp_mat_items_per_row, unsigned int sp_mat_aligned_items_per_row, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:2014
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::sliced_ell_matrix::handle3
handle_type & handle3()
Definition: sliced_ell_matrix.hpp:115
viennacl::matrix_base::row_major
bool row_major() const
Definition: matrix_def.hpp:248
viennacl::linalg::max
NumericT max(std::vector< NumericT > const &v1)
Definition: maxmin.hpp:47
viennacl::linalg::cuda::compressed_matrix_vec_mul_adaptive_kernel
__global__ void compressed_matrix_vec_mul_adaptive_kernel(const unsigned int *row_indices, const unsigned int *column_indices, const unsigned int *row_blocks, const NumericT *elements, unsigned int num_blocks, const NumericT *x, unsigned int start_x, unsigned int inc_x, NumericT alpha, NumericT *result, unsigned int start_result, unsigned int inc_result, unsigned int size_result, NumericT beta)
Definition: sparse_matrix_operations.hpp:182
viennacl::linalg::cuda::coordinate_matrix_d_tr_mat_mul_kernel
__global__ void coordinate_matrix_d_tr_mat_mul_kernel(const unsigned int *coords, const NumericT *elements, const unsigned int *group_boundaries, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:1547
viennacl::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
viennacl::compressed_compressed_matrix::nnz1
const vcl_size_t & nnz1() const
Returns the number of nonzero entries.
Definition: compressed_compressed_matrix.hpp:464
viennacl::hyb_matrix::ell_nnz
vcl_size_t ell_nnz() const
Definition: hyb_matrix.hpp:102
viennacl::linalg::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::linalg::cuda::compressed_matrix_d_mat_mul_kernel
__global__ void compressed_matrix_d_mat_mul_kernel(const unsigned int *sp_mat_row_indices, const unsigned int *sp_mat_col_indices, const NumericT *sp_mat_elements, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:430
viennacl::op_trans
A tag class representing transposed matrices.
Definition: forwards.h:220
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
VIENNACL_CUDA_LAST_ERROR_CHECK
#define VIENNACL_CUDA_LAST_ERROR_CHECK(message)
Definition: common.hpp:30
viennacl::column_major
A tag for column-major storage of a dense matrix.
Definition: forwards.h:321
viennacl::cuda_arg
NumericT * cuda_arg(scalar< NumericT > &obj)
Convenience helper function for extracting the CUDA handle from a ViennaCL scalar. Non-const version.
Definition: common.hpp:39
viennacl::linalg::cuda::detail::spmv_alpha_beta
Definition: sparse_matrix_operations.hpp:127
viennacl::hyb_matrix::handle5
const handle_type & handle5() const
Definition: hyb_matrix.hpp:109
viennacl::hyb_matrix::size1
vcl_size_t size1() const
Definition: hyb_matrix.hpp:98
viennacl::matrix_expression::lhs
LHS & lhs() const
Get left hand side operand.
Definition: matrix.hpp:66
viennacl::vector_base::start
size_type start() const
Returns the offset within the buffer.
Definition: vector_def.hpp:122
viennacl::compressed_matrix::blocks1
const vcl_size_t & blocks1() const
Returns the internal number of row blocks for an adaptive SpMV.
Definition: compressed_matrix.hpp:933
viennacl::ell_matrix::internal_maxnnz
vcl_size_t internal_maxnnz() const
Definition: ell_matrix.hpp:94
viennacl::sliced_ell_matrix::handle1
handle_type & handle1()
Definition: sliced_ell_matrix.hpp:109
scalar.hpp
Implementation of the ViennaCL scalar class.
vector_operations.hpp
Implementations of NMF operations using CUDA.
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864
viennacl::linalg::min
NumericT min(std::vector< NumericT > const &v1)
Definition: maxmin.hpp:91
viennacl::coordinate_matrix
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...
Definition: coordinate_matrix.hpp:186
viennacl::linalg::cuda::compressed_matrix_d_tr_mat_mul_kernel
__global__ void compressed_matrix_d_tr_mat_mul_kernel(const unsigned int *sp_mat_row_indices, const unsigned int *sp_mat_col_indices, const NumericT *sp_mat_elements, const NumericT *d_mat, unsigned int d_mat_row_start, unsigned int d_mat_col_start, unsigned int d_mat_row_inc, unsigned int d_mat_col_inc, unsigned int d_mat_row_size, unsigned int d_mat_col_size, unsigned int d_mat_internal_rows, unsigned int d_mat_internal_cols, NumericT *result, unsigned int result_row_start, unsigned int result_col_start, unsigned int result_row_inc, unsigned int result_col_inc, unsigned int result_row_size, unsigned int result_col_size, unsigned int result_internal_rows, unsigned int result_internal_cols)
Definition: sparse_matrix_operations.hpp:584