linearoperator_matvec.hpp

Go to the documentation of this file.

#pragma once
 
namespace monolish {
 
namespace {
template <typename T, typename VEC1, typename VEC2>
void matvec_core(const matrix::LinearOperator<T> &A, const VEC1 &x, VEC2 &y) {
  Logger &logger = Logger::get_instance();
  logger.func_in(monolish_func);
 
  // err, M = MN * N
  assert(A.get_row() == y.size());
  assert(A.get_col() == x.size());
  assert(util::is_same_device_mem_stat(A, x, y));
  assert(A.get_matvec_init_flag());
 
  const size_t xoffset = x.get_offset();
  const size_t yoffset = y.get_offset();
 
  vector<T> x_tmp(x.size(), 0);
  if (x.get_device_mem_stat())
    x_tmp.send();
  vector<T> y_tmp(y.size(), 0);
  if (y.get_device_mem_stat())
    y_tmp.send();
 
  internal::vcopy(x.size(), x.data() + xoffset, x_tmp.data(),
                  x.get_device_mem_stat());
 
  y_tmp = A.get_matvec()(x_tmp);
 
  internal::vcopy(y.size(), y_tmp.data(), y.data() + yoffset,
                  y.get_device_mem_stat());
 
  logger.func_out();
}
 
template <typename T, typename VEC1, typename VEC2>
void rmatvec_core(const matrix::LinearOperator<T> &A, const VEC1 &x, VEC2 &y) {
  Logger &logger = Logger::get_instance();
  logger.func_in(monolish_func);
 
  // err, M = MN * N
  assert(A.get_row() == x.size());
  assert(A.get_col() == y.size());
  assert(util::is_same_device_mem_stat(A, x, y));
  assert(A.get_rmatvec_init_flag());
 
  const size_t xoffset = x.get_offset();
  const size_t yoffset = y.get_offset();
 
  vector<T> x_tmp(x.size(), 0);
  if (x.get_device_mem_stat())
    x_tmp.send();
  vector<T> y_tmp(y.size(), 0);
  if (y.get_device_mem_stat())
    y_tmp.send();
 
  internal::vcopy(x.size(), x.data() + xoffset, x_tmp.data(),
                  x.get_device_mem_stat());
 
  y_tmp = A.get_rmatvec()(x_tmp);
 
  internal::vcopy(y.size(), y_tmp.data(), y.data() + yoffset,
                  y.get_device_mem_stat());
 
  logger.func_out();
}
 
} // namespace
 
} // namespace monolish