matrix-free.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>

//
// ViennaCL includes
//
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/cg.hpp"
#include "viennacl/linalg/bicgstab.hpp"
#include "viennacl/linalg/gmres.hpp"

template<typename NumericT>
class MyOperator
{
public:
  MyOperator(std::size_t N) : N_(N) {}

  // Dispatcher for y = Ax
  void apply(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const
  {
#if defined(VIENNACL_WITH_CUDA)
    if (viennacl::traits::active_handle_id(x) == viennacl::CUDA_MEMORY)
      apply_cuda(x, y);
#endif

#if defined(VIENNACL_WITH_OPENCL)
    if (viennacl::traits::active_handle_id(x) == viennacl::OPENCL_MEMORY)
      apply_opencl(x, y);
#endif

    if (viennacl::traits::active_handle_id(x) == viennacl::MAIN_MEMORY)
      apply_host(x, y);
  }

  std::size_t size1() const { return N_ * N_; }

private:

#if defined(VIENNACL_WITH_CUDA)
  void apply_cuda(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const;
#endif

#if defined(VIENNACL_WITH_OPENCL)
  void apply_opencl(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const;
#endif

  void apply_host(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const;

  std::size_t N_;
};


int main()
{
  typedef float       ScalarType;  // feel free to change to double (and change OpenCL kernel argument types accordingly)

  std::size_t N = 10;
  viennacl::vector<ScalarType> rhs = viennacl::scalar_vector<ScalarType>(N*N, ScalarType(-1));
  MyOperator<ScalarType> op(N);

  viennacl::vector<ScalarType> result = viennacl::linalg::solve(op, rhs, viennacl::linalg::cg_tag());

  std::cout.precision(3);
  std::cout << std::fixed;
  std::cout << "Result value map: " << std::endl;
  std::cout << std::endl << "^ y " << std::endl;
  for (std::size_t i=0; i<N; ++i)
  {
    std::cout << "|  ";
    for (std::size_t j=0; j<N; ++j)
      std::cout << result[i * N + j] << "  ";
    std::cout << std::endl;
  }
  std::cout << "*---------------------------------------------> x" << std::endl;

  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;

  return EXIT_SUCCESS;
}




template<typename NumericT>
void MyOperator<NumericT>::apply_host(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const
{
  NumericT const * values_x = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(x.handle());
  NumericT       * values_y = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(y.handle());

  NumericT dx = NumericT(1) / NumericT(N_ + 1);
  NumericT dy = NumericT(1) / NumericT(N_ + 1);

  // feel free to use
  //  #pragma omp parallel for
  // here
  for (std::size_t i=0; i<N_; ++i)
    for (std::size_t j=0; j<N_; ++j)
    {
      NumericT value_right  = (j < N_ - 1) ? values_x[ i   *N_ + j + 1] : 0;
      NumericT value_left   = (j > 0     ) ? values_x[ i   *N_ + j - 1] : 0;
      NumericT value_top    = (i < N_ - 1) ? values_x[(i+1)*N_ + j    ] : 0;
      NumericT value_bottom = (i > 0     ) ? values_x[(i-1)*N_ + j    ] : 0;
      NumericT value_center = values_x[i*N_ + j];

      values_y[i*N_ + j] =   ((value_right - value_center) / dx - (value_center - value_left)   / dx) / dx
                           + ((value_top   - value_center) / dy - (value_center - value_bottom) / dy) / dy;
    }
}


#if defined(VIENNACL_WITH_CUDA)
template<typename NumericT>
__global__ void apply_cuda_kernel(NumericT const * values_x,
                                  NumericT       * values_y,
                                  std::size_t N)
{
  NumericT dx = NumericT(1) / (N + 1);
  NumericT dy = NumericT(1) / (N + 1);

  for (std::size_t i = blockIdx.x; i < N; i += gridDim.x)
    for (std::size_t j = threadIdx.x; j < N; j += blockDim.x)
    {
      NumericT value_right  = (j < N - 1) ? values_x[ i   *N + j + 1] : 0;
      NumericT value_left   = (j > 0    ) ? values_x[ i   *N + j - 1] : 0;
      NumericT value_top    = (i < N - 1) ? values_x[(i+1)*N + j    ] : 0;
      NumericT value_bottom = (i > 0    ) ? values_x[(i-1)*N + j    ] : 0;
      NumericT value_center = values_x[i*N + j];

      values_y[i*N + j] =   ((value_right - value_center) / dx - (value_center - value_left)   / dx) / dx
                          + ((value_top   - value_center) / dy - (value_center - value_bottom) / dy) / dy;
    }
}
#endif


#if defined(VIENNACL_WITH_CUDA)
template<typename NumericT>
void MyOperator<NumericT>::apply_cuda(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const
{
  apply_cuda_kernel<<<128, 128>>>(viennacl::cuda_arg(x), viennacl::cuda_arg(y), N_);
}
#endif




#if defined(VIENNACL_WITH_OPENCL)
static const char * my_compute_program =
"typedef float NumericT; \n"
"__kernel void apply_opencl_kernel(__global NumericT const * values_x, \n"
"                                  __global NumericT       * values_y, \n"
"                                  unsigned int N) {\n"

"      NumericT dx = (NumericT)1 / (N + 1); \n"
"      NumericT dy = (NumericT)1 / (N + 1); \n"

"      for (unsigned int i = get_group_id(0); i < N; i += get_num_groups(0)) \n"
"        for (unsigned int j = get_local_id(0); j < N; j += get_local_size(0)) { \n"

"          NumericT value_right  = (j < N - 1) ? values_x[ i   *N + j + 1] : 0; \n"
"          NumericT value_left   = (j > 0    ) ? values_x[ i   *N + j - 1] : 0; \n"
"          NumericT value_top    = (i < N - 1) ? values_x[(i+1)*N + j    ] : 0; \n"
"          NumericT value_bottom = (i > 0    ) ? values_x[(i-1)*N + j    ] : 0; \n"
"          NumericT value_center = values_x[i*N + j]; \n"

"          values_y[i*N + j] =   ((value_right - value_center) / dx - (value_center - value_left)   / dx) / dx  \n"
"                              + ((value_top   - value_center) / dy - (value_center - value_bottom) / dy) / dy; \n"
"        }  \n"
"    } \n";
#endif

#if defined(VIENNACL_WITH_OPENCL)
template<typename NumericT>
void MyOperator<NumericT>::apply_opencl(viennacl::vector_base<NumericT> const & x, viennacl::vector_base<NumericT> & y) const
  {
    viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
    static bool first_run = true;
    if (first_run) {
      ctx.add_program(my_compute_program, "my_compute_program");
      first_run = false;
    }
    viennacl::ocl::kernel & my_kernel = ctx.get_kernel("my_compute_program", "apply_opencl_kernel");

    viennacl::ocl::enqueue(my_kernel(x, y, static_cast<cl_uint>(N_)));
  }
#endif