All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
direct_solve.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_CUDA_DIRECT_SOLVE_HPP
#define VIENNACL_LINALG_CUDA_DIRECT_SOLVE_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/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/cuda/common.hpp"
namespace viennacl
{
namespace linalg
{
namespace cuda
{
template<typename NumericT>
__global__ void matrix_matrix_upper_solve_kernel(
const NumericT * A,
unsigned int A_start1, unsigned int A_start2,
unsigned int A_inc1, unsigned int A_inc2,
unsigned int A_size1, unsigned int A_size2,
unsigned int A_internal_size1, unsigned int A_internal_size2,
bool row_major_A,
NumericT * B,
unsigned int B_start1, unsigned int B_start2,
unsigned int B_inc1, unsigned int B_inc2,
unsigned int B_size1, unsigned int B_size2,
unsigned int B_internal_size1, unsigned int B_internal_size2,
bool row_major_B,
bool unit_diagonal)
{
NumericT temp;
NumericT entry_A;
for (unsigned int row_cnt = 0; row_cnt < A_size1; ++row_cnt)
{
unsigned int row = A_size1 - 1 - row_cnt;
if (!unit_diagonal)
{
__syncthreads();
if (threadIdx.x == 0)
{
if (row_major_B)
B[(row * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)] /= (row_major_A) ? A[(row * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)]
: A[(row * A_inc1 + A_start1) + (row * A_inc2 + A_start2)*A_internal_size1];
else //if (!row_major_B)
B[(row * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1] /= (row_major_A) ? A[(row * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)]
: A[(row * A_inc1 + A_start1) + (row * A_inc2 + A_start2)*A_internal_size1];
}
}
__syncthreads();
if (row_major_B)
temp = B[(row * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)];
else //if (!row_major_B)
temp = B[(row * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1];
//eliminate column of op(A) with index 'row' in parallel: " << std::endl;
for (unsigned int elim = threadIdx.x; elim < row; elim += blockDim.x)
{
if (row_major_A)
entry_A = A[(elim * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)];
else //if (!row_major_A)
entry_A = A[(elim * A_inc1 + A_start1) + (row * A_inc2 + A_start2) * A_internal_size1];
if (row_major_B)
B[(elim * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)] -= temp * entry_A;
else //if (!row_major_B)
B[(elim * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1] -= temp * entry_A;
}
}
}
template<typename NumericT>
__global__ void matrix_matrix_lower_solve_kernel(
const NumericT * A,
unsigned int A_start1, unsigned int A_start2,
unsigned int A_inc1, unsigned int A_inc2,
unsigned int A_size1, unsigned int A_size2,
unsigned int A_internal_size1, unsigned int A_internal_size2,
bool row_major_A,
NumericT * B,
unsigned int B_start1, unsigned int B_start2,
unsigned int B_inc1, unsigned int B_inc2,
unsigned int B_size1, unsigned int B_size2,
unsigned int B_internal_size1, unsigned int B_internal_size2,
bool row_major_B,
bool unit_diagonal)
{
NumericT temp;
NumericT entry_A;
for (unsigned int row = 0; row < A_size1; ++row)
{
if (!unit_diagonal)
{
__syncthreads();
if (threadIdx.x == 0)
{
if (row_major_B)
B[(row * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)] /= (row_major_A) ? A[(row * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)]
: A[(row * A_inc1 + A_start1) + (row * A_inc2 + A_start2)*A_internal_size1];
else //if (!row_major_B)
B[(row * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1] /= (row_major_A) ? A[(row * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)]
: A[(row * A_inc1 + A_start1) + (row * A_inc2 + A_start2)*A_internal_size1];
}
}
__syncthreads();
if (row_major_B)
temp = B[(row * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)];
else //if (!row_major_B)
temp = B[(row * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1];
//eliminate column of op(A) with index 'row' in parallel: " << std::endl;
for (unsigned int elim = row + threadIdx.x + 1; elim < A_size1; elim += blockDim.x)
{
if (row_major_A)
entry_A = A[(elim * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)];
else //if (!row_major_A)
entry_A = A[(elim * A_inc1 + A_start1) + (row * A_inc2 + A_start2) * A_internal_size1];
if (row_major_B)
B[(elim * B_inc1 + B_start1) * B_internal_size2 + (blockIdx.x * B_inc2 + B_start2)] -= temp * entry_A;
else //if (!row_major_B)
B[(elim * B_inc1 + B_start1) + (blockIdx.x * B_inc2 + B_start2) * B_internal_size1] -= temp * entry_A;
}
}
}
namespace detail
{
template<typename TagT>
bool is_unit_solve(TagT const & tag) { return false; }
inline bool is_unit_solve(viennacl::linalg::unit_lower_tag) { return true; }
inline bool is_unit_solve(viennacl::linalg::unit_upper_tag) { return true; }
template<typename TagT>
bool is_upper_solve(TagT const & tag) { return false; }
inline bool is_upper_solve(viennacl::linalg::upper_tag) { return true; }
inline bool is_upper_solve(viennacl::linalg::unit_upper_tag) { return true; }
template<typename Matrix1T, typename Matrix2T, typename SolverTagT>
void inplace_solve_impl(Matrix1T const & A,
Matrix2T & B,
SolverTagT const & tag)
{
typedef typename viennacl::result_of::cpu_value_type<Matrix1T>::type value_type;
dim3 threads(128);
dim3 grid(B.size2());
if (is_upper_solve(tag))
{
matrix_matrix_upper_solve_kernel<<<grid,threads>>>(viennacl::cuda_arg(A),
static_cast<unsigned int>(viennacl::traits::start1(A)), static_cast<unsigned int>(viennacl::traits::start2(A)),
static_cast<unsigned int>(viennacl::traits::stride1(A)), static_cast<unsigned int>(viennacl::traits::stride2(A)),
static_cast<unsigned int>(viennacl::traits::size1(A)), static_cast<unsigned int>(viennacl::traits::size2(A)),
static_cast<unsigned int>(viennacl::traits::internal_size1(A)), static_cast<unsigned int>(viennacl::traits::internal_size2(A)),
bool(A.row_major()),
viennacl::cuda_arg(B),
static_cast<unsigned int>(viennacl::traits::start1(B)), static_cast<unsigned int>(viennacl::traits::start2(B)),
static_cast<unsigned int>(viennacl::traits::stride1(B)), static_cast<unsigned int>(viennacl::traits::stride2(B)),
static_cast<unsigned int>(viennacl::traits::size1(B)), static_cast<unsigned int>(viennacl::traits::size2(B)),
static_cast<unsigned int>(viennacl::traits::internal_size1(B)), static_cast<unsigned int>(viennacl::traits::internal_size2(B)),
bool(B.row_major()),
is_unit_solve(tag)
);
}
else
{
matrix_matrix_lower_solve_kernel<<<grid,threads>>>(viennacl::cuda_arg(A),
static_cast<unsigned int>(viennacl::traits::start1(A)), static_cast<unsigned int>(viennacl::traits::start2(A)),
static_cast<unsigned int>(viennacl::traits::stride1(A)), static_cast<unsigned int>(viennacl::traits::stride2(A)),
static_cast<unsigned int>(viennacl::traits::size1(A)), static_cast<unsigned int>(viennacl::traits::size2(A)),
static_cast<unsigned int>(viennacl::traits::internal_size1(A)), static_cast<unsigned int>(viennacl::traits::internal_size2(A)),
bool(A.row_major()),
viennacl::cuda_arg(B),
static_cast<unsigned int>(viennacl::traits::start1(B)), static_cast<unsigned int>(viennacl::traits::start2(B)),
static_cast<unsigned int>(viennacl::traits::stride1(B)), static_cast<unsigned int>(viennacl::traits::stride2(B)),
static_cast<unsigned int>(viennacl::traits::size1(B)), static_cast<unsigned int>(viennacl::traits::size2(B)),
static_cast<unsigned int>(viennacl::traits::internal_size1(B)), static_cast<unsigned int>(viennacl::traits::internal_size2(B)),
bool(B.row_major()),
is_unit_solve(tag)
);
}
}
}
//
// Note: By convention, all size checks are performed in the calling frontend. No need to double-check here.
//
template<typename NumericT, typename SolverTagT>
void inplace_solve(matrix_base<NumericT> const & A,
matrix_base<NumericT> & B,
SolverTagT tag)
{
detail::inplace_solve_impl(A, B, tag);
}
//
// Solve on vector
//
template<typename NumericT>
__global__ void triangular_substitute_inplace_row_kernel(
NumericT const * A,
unsigned int A_start1, unsigned int A_start2,
unsigned int A_inc1, unsigned int A_inc2,
unsigned int A_size1, unsigned int A_size2,
unsigned int A_internal_size1, unsigned int A_internal_size2,
NumericT * v,
unsigned int v_start,
unsigned int v_inc,
unsigned int v_size,
unsigned int options)
{
NumericT temp;
unsigned int unit_diagonal_flag = (options & (1 << 0));
unsigned int is_lower_solve = (options & (1 << 2));
unsigned int row;
for (unsigned int rows_processed = 0; rows_processed < A_size1; ++rows_processed) //Note: A required to be square
{
row = is_lower_solve ? rows_processed : ((A_size1 - rows_processed) - 1);
if (!unit_diagonal_flag)
{
__syncthreads();
if (threadIdx.x == 0)
v[row * v_inc + v_start] /= A[(row * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)];
}
__syncthreads();
temp = v[row * v_inc + v_start];
for (int elim = (is_lower_solve ? (row + threadIdx.x + 1) : threadIdx.x);
elim < (is_lower_solve ? A_size1 : row);
elim += blockDim.x)
v[elim * v_inc + v_start] -= temp * A[(elim * A_inc1 + A_start1) * A_internal_size2 + (row * A_inc2 + A_start2)];
}
}
template<typename NumericT>
__global__ void triangular_substitute_inplace_col_kernel(
NumericT const * A,
unsigned int A_start1, unsigned int A_start2,
unsigned int A_inc1, unsigned int A_inc2,
unsigned int A_size1, unsigned int A_size2,
unsigned int A_internal_size1, unsigned int A_internal_size2,
NumericT * v,
unsigned int v_start,
unsigned int v_inc,
unsigned int v_size,
unsigned int options)
{
NumericT temp;
unsigned int unit_diagonal_flag = (options & (1 << 0));
unsigned int is_lower_solve = (options & (1 << 2));
unsigned int row;
for (unsigned int rows_processed = 0; rows_processed < A_size1; ++rows_processed) //Note: A required to be square
{
row = is_lower_solve ? rows_processed : ((A_size1 - rows_processed) - 1);
if (!unit_diagonal_flag)
{
__syncthreads();
if (threadIdx.x == 0)
v[row * v_inc + v_start] /= A[(row * A_inc1 + A_start1) + (row * A_inc2 + A_start2) * A_internal_size1];
}
__syncthreads();
temp = v[row * v_inc + v_start];
for (int elim = (is_lower_solve ? (row + threadIdx.x + 1) : threadIdx.x);
elim < (is_lower_solve ? A_size1 : row);
elim += blockDim.x)
v[elim * v_inc + v_start] -= temp * A[(elim * A_inc1 + A_start1) + (row * A_inc2 + A_start2) * A_internal_size1];
}
}
namespace detail
{
inline unsigned int get_option_for_solver_tag(viennacl::linalg::upper_tag) { return 0; }
inline unsigned int get_option_for_solver_tag(viennacl::linalg::unit_upper_tag) { return (1 << 0); }
inline unsigned int get_option_for_solver_tag(viennacl::linalg::lower_tag) { return (1 << 2); }
inline unsigned int get_option_for_solver_tag(viennacl::linalg::unit_lower_tag) { return (1 << 2) | (1 << 0); }
template<typename MatrixT, typename VectorT>
void inplace_solve_vector_impl(MatrixT const & mat,
VectorT & vec,
unsigned int options)
{
typedef typename viennacl::result_of::cpu_value_type<MatrixT>::type value_type;
if (mat.row_major())
{
triangular_substitute_inplace_row_kernel<<<1, 128>>>(viennacl::cuda_arg(mat),
static_cast<unsigned int>(viennacl::traits::start1(mat)), static_cast<unsigned int>(viennacl::traits::start2(mat)),
static_cast<unsigned int>(viennacl::traits::stride1(mat)), static_cast<unsigned int>(viennacl::traits::stride2(mat)),
static_cast<unsigned int>(viennacl::traits::size1(mat)), static_cast<unsigned int>(viennacl::traits::size2(mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(mat)),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(viennacl::traits::start(vec)),
static_cast<unsigned int>(viennacl::traits::stride(vec)),
static_cast<unsigned int>(viennacl::traits::size(vec)),
options
);
}
else
{
triangular_substitute_inplace_col_kernel<<<1, 128>>>(viennacl::cuda_arg(mat),
static_cast<unsigned int>(viennacl::traits::start1(mat)), static_cast<unsigned int>(viennacl::traits::start2(mat)),
static_cast<unsigned int>(viennacl::traits::stride1(mat)), static_cast<unsigned int>(viennacl::traits::stride2(mat)),
static_cast<unsigned int>(viennacl::traits::size1(mat)), static_cast<unsigned int>(viennacl::traits::size2(mat)),
static_cast<unsigned int>(viennacl::traits::internal_size1(mat)), static_cast<unsigned int>(viennacl::traits::internal_size2(mat)),
viennacl::cuda_arg(vec),
static_cast<unsigned int>(viennacl::traits::start(vec)),
static_cast<unsigned int>(viennacl::traits::stride(vec)),
static_cast<unsigned int>(viennacl::traits::size(vec)),
options
);
}
}
}
template<typename NumericT, typename SolverTagT>
void inplace_solve(matrix_base<NumericT> const & mat,
vector_base<NumericT> & vec,
SolverTagT)
{
unsigned int options = detail::get_option_for_solver_tag(SolverTagT());
detail::inplace_solve_vector_impl(mat, vec, options);
}
}
}
}
#endif
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::linalg::cuda::detail::inplace_solve_vector_impl
void inplace_solve_vector_impl(MatrixT const &mat, VectorT &vec, unsigned int options)
Definition: direct_solve.hpp:354
viennacl::traits::stride1
result_of::size_type< matrix_base< NumericT > >::type stride1(matrix_base< NumericT > const &s)
Definition: stride.hpp:55
viennacl::linalg::cuda::triangular_substitute_inplace_row_kernel
__global__ void triangular_substitute_inplace_row_kernel(NumericT const *A, unsigned int A_start1, unsigned int A_start2, unsigned int A_inc1, unsigned int A_inc2, unsigned int A_size1, unsigned int A_size2, unsigned int A_internal_size1, unsigned int A_internal_size2, NumericT *v, unsigned int v_start, unsigned int v_inc, unsigned int v_size, unsigned int options)
Definition: direct_solve.hpp:266
matrix.hpp
Implementation of the dense matrix class.
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::linalg::cuda::matrix_matrix_lower_solve_kernel
__global__ void matrix_matrix_lower_solve_kernel(const NumericT *A, unsigned int A_start1, unsigned int A_start2, unsigned int A_inc1, unsigned int A_inc2, unsigned int A_size1, unsigned int A_size2, unsigned int A_internal_size1, unsigned int A_internal_size2, bool row_major_A, NumericT *B, unsigned int B_start1, unsigned int B_start2, unsigned int B_inc1, unsigned int B_inc2, unsigned int B_size1, unsigned int B_size2, unsigned int B_internal_size1, unsigned int B_internal_size2, bool row_major_B, bool unit_diagonal)
Definition: direct_solve.hpp:107
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::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::matrix_base< NumericT >
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::traits::start1
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:65
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::linalg::cuda::matrix_matrix_upper_solve_kernel
__global__ void matrix_matrix_upper_solve_kernel(const NumericT *A, unsigned int A_start1, unsigned int A_start2, unsigned int A_inc1, unsigned int A_inc2, unsigned int A_size1, unsigned int A_size2, unsigned int A_internal_size1, unsigned int A_internal_size2, bool row_major_A, NumericT *B, unsigned int B_start1, unsigned int B_start2, unsigned int B_inc1, unsigned int B_inc2, unsigned int B_size1, unsigned int B_size2, unsigned int B_internal_size1, unsigned int B_internal_size2, bool row_major_B, bool unit_diagonal)
Definition: direct_solve.hpp:41
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::cuda::detail::is_unit_solve
bool is_unit_solve(TagT const &tag)
Definition: direct_solve.hpp:177
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::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::traits::start
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:44
viennacl::linalg::cuda::triangular_substitute_inplace_col_kernel
__global__ void triangular_substitute_inplace_col_kernel(NumericT const *A, unsigned int A_start1, unsigned int A_start2, unsigned int A_inc1, unsigned int A_inc2, unsigned int A_size1, unsigned int A_size2, unsigned int A_internal_size1, unsigned int A_internal_size2, NumericT *v, unsigned int v_start, unsigned int v_inc, unsigned int v_size, unsigned int options)
Definition: direct_solve.hpp:307
viennacl::vector_base< NumericT >
viennacl::result_of::cpu_value_type::type
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
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::linalg::cuda::detail::inplace_solve_impl
void inplace_solve_impl(Matrix1T const &A, Matrix2T &B, SolverTagT const &tag)
Definition: direct_solve.hpp:189
common.hpp
Common routines for CUDA execution.
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::linalg::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
viennacl::linalg::cuda::detail::get_option_for_solver_tag
unsigned int get_option_for_solver_tag(viennacl::linalg::upper_tag)
Definition: direct_solve.hpp:348
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::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864
viennacl::linalg::cuda::detail::is_upper_solve
bool is_upper_solve(TagT const &tag)
Definition: direct_solve.hpp:183