All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
libviennacl_blas2.cpp
Go to the documentation of this file.
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the PDF manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
// include necessary system headers
#include <iostream>
#include <vector>
// Some helper functions for this tutorial:
#include "viennacl.hpp"
#include "viennacl/vector.hpp"
template<typename ScalarType>
ScalarType diff(ScalarType const & s1, ScalarType const & s2)
{
if (s1 > s2 || s1 < s2)
return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
return ScalarType(0);
}
template<typename ScalarType, typename ViennaCLVectorType>
ScalarType diff(std::vector<ScalarType> const & v1, ViennaCLVectorType const & vcl_vec)
{
std::vector<ScalarType> v2_cpu(vcl_vec.size());
viennacl::backend::finish();
viennacl::copy(vcl_vec, v2_cpu);
ScalarType inf_norm = 0;
for (unsigned int i=0;i<v1.size(); ++i)
{
if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
else
v2_cpu[i] = 0.0;
if (v2_cpu[i] > inf_norm)
inf_norm = v2_cpu[i];
}
return inf_norm;
}
template<typename T, typename U, typename EpsilonT>
void check(T const & t, U const & u, EpsilonT eps)
{
EpsilonT rel_error = std::fabs(static_cast<EpsilonT>(diff(t,u)));
if (rel_error > eps)
{
std::cerr << "Relative error: " << rel_error << std::endl;
std::cerr << "Aborting!" << std::endl;
exit(EXIT_FAILURE);
}
std::cout << "SUCCESS ";
}
int main()
{
std::size_t size1 = 13; // at least 7
std::size_t size2 = 11; // at least 7
float eps_float = 1e-5f;
double eps_double = 1e-12;
ViennaCLBackend my_backend;
ViennaCLBackendCreate(&my_backend);
std::vector<float> ref_float_x(size1); for (std::size_t i=0; i<size1; ++i) ref_float_x[i] = static_cast<float>(i);
std::vector<float> ref_float_y(size2); for (std::size_t i=0; i<size2; ++i) ref_float_y[i] = static_cast<float>(size2 - i);
std::vector<float> ref_float_A(size1*size2); for (std::size_t i=0; i<size1*size2; ++i) ref_float_A[i] = static_cast<float>(3*i);
std::vector<float> ref_float_B(size1*size2); for (std::size_t i=0; i<size1*size2; ++i) ref_float_B[i] = static_cast<float>(2*i);
std::vector<double> ref_double_x(size1, 1.0); for (std::size_t i=0; i<size1; ++i) ref_double_x[i] = static_cast<double>(i);
std::vector<double> ref_double_y(size2, 2.0); for (std::size_t i=0; i<size2; ++i) ref_double_y[i] = static_cast<double>(size2 - i);
std::vector<double> ref_double_A(size1*size2, 3.0); for (std::size_t i=0; i<size1*size2; ++i) ref_double_A[i] = static_cast<double>(3*i);
std::vector<double> ref_double_B(size1*size2, 4.0); for (std::size_t i=0; i<size1*size2; ++i) ref_double_B[i] = static_cast<double>(2*i);
// Host setup
viennacl::vector<float> host_float_x = viennacl::scalar_vector<float>(size1, 1.0f, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1; ++i) host_float_x[i] = float(i);
viennacl::vector<float> host_float_y = viennacl::scalar_vector<float>(size2, 2.0f, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size2; ++i) host_float_y[i] = float(size2 - i);
viennacl::vector<float> host_float_A = viennacl::scalar_vector<float>(size1*size2, 3.0f, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) host_float_A[i] = float(3*i);
viennacl::vector<float> host_float_B = viennacl::scalar_vector<float>(size1*size2, 4.0f, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) host_float_B[i] = float(2*i);
viennacl::vector<double> host_double_x = viennacl::scalar_vector<double>(size1, 1.0, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1; ++i) host_double_x[i] = double(i);
viennacl::vector<double> host_double_y = viennacl::scalar_vector<double>(size2, 2.0, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size2; ++i) host_double_y[i] = double(size2 - i);
viennacl::vector<double> host_double_A = viennacl::scalar_vector<double>(size1*size2, 3.0, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) host_double_A[i] = double(3*i);
viennacl::vector<double> host_double_B = viennacl::scalar_vector<double>(size1*size2, 4.0, viennacl::context(viennacl::MAIN_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) host_double_B[i] = double(2*i);
// CUDA setup
#ifdef VIENNACL_WITH_CUDA
viennacl::vector<float> cuda_float_x = viennacl::scalar_vector<float>(size1, 1.0f, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1; ++i) cuda_float_x[i] = float(i);
viennacl::vector<float> cuda_float_y = viennacl::scalar_vector<float>(size2, 2.0f, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size2; ++i) cuda_float_y[i] = float(size2 - i);
viennacl::vector<float> cuda_float_A = viennacl::scalar_vector<float>(size1*size2, 3.0f, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) cuda_float_A[i] = float(3*i);
viennacl::vector<float> cuda_float_B = viennacl::scalar_vector<float>(size1*size2, 4.0f, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) cuda_float_B[i] = float(2*i);
viennacl::vector<double> cuda_double_x = viennacl::scalar_vector<double>(size1, 1.0, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1; ++i) cuda_double_x[i] = double(i);
viennacl::vector<double> cuda_double_y = viennacl::scalar_vector<double>(size2, 2.0, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size2; ++i) cuda_double_y[i] = double(size2 - i);
viennacl::vector<double> cuda_double_A = viennacl::scalar_vector<double>(size1*size2, 3.0, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) cuda_double_A[i] = double(3*i);
viennacl::vector<double> cuda_double_B = viennacl::scalar_vector<double>(size1*size2, 4.0, viennacl::context(viennacl::CUDA_MEMORY)); for (std::size_t i=0; i<size1*size2; ++i) cuda_double_B[i] = double(2*i);
#endif
// OpenCL setup
#ifdef VIENNACL_WITH_OPENCL
ViennaCLInt context_id = 0;
viennacl::vector<float> opencl_float_x = viennacl::scalar_vector<float>(size1, 1.0f, viennacl::context(viennacl::ocl::get_context(context_id))); for (std::size_t i=0; i<size1; ++i) opencl_float_x[i] = float(i);
viennacl::vector<float> opencl_float_y = viennacl::scalar_vector<float>(size2, 2.0f, viennacl::context(viennacl::ocl::get_context(context_id))); for (std::size_t i=0; i<size2; ++i) opencl_float_y[i] = float(size2 - i);
viennacl::vector<float> opencl_float_A = viennacl::scalar_vector<float>(size1*size2, 3.0f, viennacl::context(viennacl::ocl::get_context(context_id))); for (std::size_t i=0; i<size1*size2; ++i) opencl_float_A[i] = float(3*i);
viennacl::vector<float> opencl_float_B = viennacl::scalar_vector<float>(size1*size2, 4.0f, viennacl::context(viennacl::ocl::get_context(context_id))); for (std::size_t i=0; i<size1*size2; ++i) opencl_float_B[i] = float(2*i);
viennacl::vector<double> *opencl_double_x = NULL;
viennacl::vector<double> *opencl_double_y = NULL;
viennacl::vector<double> *opencl_double_A = NULL;
viennacl::vector<double> *opencl_double_B = NULL;
if ( viennacl::ocl::current_device().double_support() )
{
opencl_double_x = new viennacl::vector<double>(viennacl::scalar_vector<double>(size1, 1.0, viennacl::context(viennacl::ocl::get_context(context_id)))); for (std::size_t i=0; i<size1; ++i) (*opencl_double_x)[i] = double(i);
opencl_double_y = new viennacl::vector<double>(viennacl::scalar_vector<double>(size2, 2.0, viennacl::context(viennacl::ocl::get_context(context_id)))); for (std::size_t i=0; i<size2; ++i) (*opencl_double_y)[i] = double(size2 - i);
opencl_double_A = new viennacl::vector<double>(viennacl::scalar_vector<double>(size1*size2, 3.0, viennacl::context(viennacl::ocl::get_context(context_id)))); for (std::size_t i=0; i<size1*size2; ++i) (*opencl_double_A)[i] = double(3*i);
opencl_double_B = new viennacl::vector<double>(viennacl::scalar_vector<double>(size1*size2, 4.0, viennacl::context(viennacl::ocl::get_context(context_id)))); for (std::size_t i=0; i<size1*size2; ++i) (*opencl_double_B)[i] = double(2*i);
}
ViennaCLBackendSetOpenCLContextID(my_backend, context_id);
#endif
// consistency checks:
check(ref_float_x, host_float_x, eps_float);
check(ref_float_y, host_float_y, eps_float);
check(ref_float_A, host_float_A, eps_float);
check(ref_float_B, host_float_B, eps_float);
check(ref_double_x, host_double_x, eps_double);
check(ref_double_y, host_double_y, eps_double);
check(ref_double_A, host_double_A, eps_double);
check(ref_double_B, host_double_B, eps_double);
#ifdef VIENNACL_WITH_CUDA
check(ref_float_x, cuda_float_x, eps_float);
check(ref_float_y, cuda_float_y, eps_float);
check(ref_float_A, cuda_float_A, eps_float);
check(ref_float_B, cuda_float_B, eps_float);
check(ref_double_x, cuda_double_x, eps_double);
check(ref_double_y, cuda_double_y, eps_double);
check(ref_double_A, cuda_double_A, eps_double);
check(ref_double_B, cuda_double_B, eps_double);
#endif
#ifdef VIENNACL_WITH_OPENCL
check(ref_float_x, opencl_float_x, eps_float);
check(ref_float_y, opencl_float_y, eps_float);
check(ref_float_A, opencl_float_A, eps_float);
check(ref_float_B, opencl_float_B, eps_float);
if ( viennacl::ocl::current_device().double_support() )
{
check(ref_double_x, *opencl_double_x, eps_double);
check(ref_double_y, *opencl_double_y, eps_double);
check(ref_double_A, *opencl_double_A, eps_double);
check(ref_double_B, *opencl_double_B, eps_double);
}
#endif
// GEMV
std::cout << std::endl << "-- Testing xGEMV...";
for (std::size_t i=0; i<size1/3; ++i)
{
ref_float_x[i * 2 + 1] *= 0.1234f;
ref_double_x[i * 2 + 1] *= 0.1234;
for (std::size_t j=0; j<size2/4; ++j)
{
ref_float_x[i * 2 + 1] += 3.1415f * ref_float_A[(2*i+2) * size2 + 3 * j + 1] * ref_float_y[j * 3 + 1];
ref_double_x[i * 2 + 1] += 3.1415 * ref_double_A[(2*i+2) * size2 + 3 * j + 1] * ref_double_y[j * 3 + 1];
}
}
std::cout << std::endl << "Host: ";
ViennaCLHostSgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415f, viennacl::linalg::host_based::detail::extract_raw_pointer<float>(host_float_A), 2, 1, 2, 3, ViennaCLInt(size2),
viennacl::linalg::host_based::detail::extract_raw_pointer<float>(host_float_y), 1, 3,
0.1234f,
viennacl::linalg::host_based::detail::extract_raw_pointer<float>(host_float_x), 1, 2);
check(ref_float_x, host_float_x, eps_float);
ViennaCLHostDgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415, viennacl::linalg::host_based::detail::extract_raw_pointer<double>(host_double_A), 2, 1, 2, 3, ViennaCLInt(size2),
viennacl::linalg::host_based::detail::extract_raw_pointer<double>(host_double_y), 1, 3,
0.1234,
viennacl::linalg::host_based::detail::extract_raw_pointer<double>(host_double_x), 1, 2);
check(ref_double_x, host_double_x, eps_double);
#ifdef VIENNACL_WITH_CUDA
std::cout << std::endl << "CUDA: ";
ViennaCLCUDASgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415f, viennacl::cuda_arg(cuda_float_A), 2, 1, 2, 3, size2,
viennacl::cuda_arg(cuda_float_y), 1, 3,
0.1234f,
viennacl::cuda_arg(cuda_float_x), 1, 2);
check(ref_float_x, cuda_float_x, eps_float);
ViennaCLCUDADgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415, viennacl::cuda_arg(cuda_double_A), 2, 1, 2, 3, size2,
viennacl::cuda_arg(cuda_double_y), 1, 3,
0.1234,
viennacl::cuda_arg(cuda_double_x), 1, 2);
check(ref_double_x, cuda_double_x, eps_double);
#endif
#ifdef VIENNACL_WITH_OPENCL
std::cout << std::endl << "OpenCL: ";
ViennaCLOpenCLSgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415f, viennacl::traits::opencl_handle(opencl_float_A), 2, 1, 2, 3, ViennaCLInt(size2),
viennacl::traits::opencl_handle(opencl_float_y), 1, 3,
0.1234f,
viennacl::traits::opencl_handle(opencl_float_x), 1, 2);
check(ref_float_x, opencl_float_x, eps_float);
if ( viennacl::ocl::current_device().double_support() )
{
ViennaCLOpenCLDgemv(my_backend,
ViennaCLRowMajor, ViennaCLNoTrans,
ViennaCLInt(size1/3), ViennaCLInt(size2/4), 3.1415, viennacl::traits::opencl_handle(*opencl_double_A), 2, 1, 2, 3, ViennaCLInt(size2),
viennacl::traits::opencl_handle(*opencl_double_y), 1, 3,
0.1234,
viennacl::traits::opencl_handle(*opencl_double_x), 1, 2);
check(ref_double_x, *opencl_double_x, eps_double);
}
#endif
#ifdef VIENNACL_WITH_OPENCL
delete opencl_double_x;
delete opencl_double_y;
delete opencl_double_A;
delete opencl_double_B;
#endif
ViennaCLBackendDestroy(&my_backend);
//
// That's it.
//
std::cout << std::endl << "!!!! TEST COMPLETED SUCCESSFULLY !!!!" << std::endl;
return EXIT_SUCCESS;
}
ViennaCLBackendCreate
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLBackendCreate(ViennaCLBackend *backend)
Definition: backend.cpp:25
ViennaCLBackend_impl
Generic backend for CUDA, OpenCL, host-based stuff.
Definition: viennacl_private.hpp:54
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
ViennaCLOpenCLDgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLOpenCLDgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, double alpha, cl_mem A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, cl_mem x, ViennaCLInt offx, ViennaCLInt incx, double beta, cl_mem y, ViennaCLInt offy, ViennaCLInt incy)
Definition: blas2_opencl.cpp:60
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
ViennaCLBackendSetOpenCLContextID
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLBackendSetOpenCLContextID(ViennaCLBackend backend, ViennaCLInt context_id)
Definition: backend.cpp:32
ViennaCLHostSgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLHostSgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, float alpha, float *A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, float *x, ViennaCLInt offx, ViennaCLInt incx, float beta, float *y, ViennaCLInt offy, ViennaCLInt incy)
Definition: blas2_host.cpp:36
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
s1
viennacl::scalar< float > s1
Definition: global_variables.cpp:57
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::current_device
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
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
ViennaCLCUDASgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLCUDASgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, float alpha, float *A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, float *x, ViennaCLInt offx, ViennaCLInt incx, float beta, float *y, ViennaCLInt offy, ViennaCLInt incy)
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
ViennaCLHostDgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLHostDgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, double alpha, double *A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, double *x, ViennaCLInt offx, ViennaCLInt incx, double beta, double *y, ViennaCLInt offy, ViennaCLInt incy)
Definition: blas2_host.cpp:60
viennacl::CUDA_MEMORY
Definition: forwards.h:350
ViennaCLBackendDestroy
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLBackendDestroy(ViennaCLBackend *backend)
Definition: backend.cpp:39
main
int main()
Definition: libviennacl_blas2.cpp:75
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
ViennaCLInt
int ViennaCLInt
Definition: viennacl.hpp:48
ViennaCLCUDADgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLCUDADgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, double alpha, double *A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, double *x, ViennaCLInt offx, ViennaCLInt incx, double beta, double *y, ViennaCLInt offy, ViennaCLInt incy)
viennacl::vector
Definition: forwards.h:266
viennacl::MAIN_MEMORY
Definition: forwards.h:348
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
ViennaCLOpenCLSgemv
VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLOpenCLSgemv(ViennaCLBackend backend, ViennaCLOrder order, ViennaCLTranspose transA, ViennaCLInt m, ViennaCLInt n, float alpha, cl_mem A, ViennaCLInt offA_row, ViennaCLInt offA_col, ViennaCLInt incA_row, ViennaCLInt incA_col, ViennaCLInt lda, cl_mem x, ViennaCLInt offx, ViennaCLInt incx, float beta, cl_mem y, ViennaCLInt offy, ViennaCLInt incy)
Definition: blas2_opencl.cpp:36
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::scalar_vector
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: vector_def.hpp:87
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
ViennaCLRowMajor
Definition: viennacl.hpp:64
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
ViennaCLNoTrans
Definition: viennacl.hpp:71
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
check
void check(T const &t, U const &u, EpsilonT eps)
Definition: libviennacl_blas2.cpp:63
viennacl::ocl::get_context
viennacl::ocl::context & get_context(long i)
Convenience function for returning the current context.
Definition: backend.hpp:225
diff
ScalarType diff(ScalarType const &s1, ScalarType const &s2)
Definition: libviennacl_blas2.cpp:33
viennacl.hpp