C++ viennacl::vector_base类代码示例

    typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
    prod_impl(const SparseMatrixType & mat,
              const viennacl::vector_base<ScalarType> & vec,
                    viennacl::vector_base<ScalarType> & result)
    {
      assert( (mat.size1() == result.size()) && bool("Size check failed for compressed matrix-vector product: size1(mat) != size(result)"));
      assert( (mat.size2() == vec.size())    && bool("Size check failed for compressed matrix-vector product: size2(mat) != size(x)"));

      switch (viennacl::traits::handle(mat).get_active_handle_id())
      {
        case viennacl::MAIN_MEMORY:
          viennacl::linalg::host_based::prod_impl(mat, vec, result);
          break;
#ifdef VIENNACL_WITH_OPENCL
        case viennacl::OPENCL_MEMORY:
          viennacl::linalg::opencl::prod_impl(mat, vec, result);
          break;
#endif
#ifdef VIENNACL_WITH_CUDA
        case viennacl::CUDA_MEMORY:
          viennacl::linalg::cuda::prod_impl(mat, vec, result);
          break;
#endif
        case viennacl::MEMORY_NOT_INITIALIZED:
          throw memory_exception("not initialised!");
        default:
          throw memory_exception("not implemented");
      }
    }

      typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
      block_inplace_solve(const matrix_expression<const SparseMatrixType, const SparseMatrixType, op_trans> & mat,
                          viennacl::backend::mem_handle const & block_index_array, vcl_size_t num_blocks,
                          viennacl::vector_base<ScalarType> const & mat_diagonal,
                          viennacl::vector_base<ScalarType> & vec,
                          SOLVERTAG tag)
      {
        assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on transposed compressed matrix: size1(mat) != size2(mat)"));
        assert( (mat.size1() == vec.size())  && bool("Size check failed for transposed compressed matrix triangular solve: size1(mat) != size(x)"));

        switch (viennacl::traits::handle(mat.lhs()).get_active_handle_id())
        {
          case viennacl::MAIN_MEMORY:
            viennacl::linalg::host_based::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
            break;
  #ifdef VIENNACL_WITH_OPENCL
          case viennacl::OPENCL_MEMORY:
            viennacl::linalg::opencl::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
            break;
  #endif
  #ifdef VIENNACL_WITH_CUDA
          case viennacl::CUDA_MEMORY:
            viennacl::linalg::cuda::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
            break;
  #endif
          case viennacl::MEMORY_NOT_INITIALIZED:
            throw memory_exception("not initialised!");
          default:
            throw memory_exception("not implemented");
        }
      }

 typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
                               viennacl::vector<SCALARTYPE> >::type
 operator-(viennacl::vector_base<SCALARTYPE> & result,
           const viennacl::vector_expression< const SparseMatrixType, const viennacl::vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
 {
   assert(proxy.lhs().size1() == result.size() && bool("Dimensions for addition of sparse matrix-vector product to vector don't match!"));
   vector<SCALARTYPE> temp(proxy.lhs().size1());
   viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), temp);
   result += temp;
   return result;
 }

    void prod_impl(const viennacl::toeplitz_matrix<SCALARTYPE, ALIGNMENT> & mat,
                   const viennacl::vector_base<SCALARTYPE> & vec,
                         viennacl::vector_base<SCALARTYPE> & result)
    {
      assert(mat.size1() == result.size());
      assert(mat.size2() == vec.size());

      viennacl::vector<SCALARTYPE> tmp(vec.size() * 4); tmp.clear();
      viennacl::vector<SCALARTYPE> tmp2(vec.size() * 4);

      viennacl::vector<SCALARTYPE> tep(mat.elements().size() * 2);
      viennacl::linalg::real_to_complex(mat.elements(), tep, mat.elements().size());



      viennacl::copy(vec.begin(), vec.end(), tmp.begin());
      viennacl::linalg::real_to_complex(tmp, tmp2, vec.size() * 2);
      viennacl::linalg::convolve(tep, tmp2, tmp);
      viennacl::linalg::complex_to_real(tmp, tmp2, vec.size() * 2);
      viennacl::copy(tmp2.begin(), tmp2.begin() + static_cast<vcl_ptrdiff_t>(vec.size()), result.begin());
    }

void prod_impl(viennacl::circulant_matrix<NumericT, AlignmentV> const & mat,
               viennacl::vector_base<NumericT> const & vec,
               viennacl::vector_base<NumericT>       & result)
{
  assert(mat.size1() == result.size() && bool("Dimension mismatch"));
  assert(mat.size2() == vec.size() && bool("Dimension mismatch"));
  //result.clear();

  //Rcpp::Rcout << "prod(circulant_matrix" << ALIGNMENT << ", vector) called with internal_nnz=" << mat.internal_nnz() << std::endl;

  viennacl::vector<NumericT> circ(mat.elements().size() * 2);
  viennacl::linalg::real_to_complex(mat.elements(), circ, mat.elements().size());

  viennacl::vector<NumericT> tmp(vec.size() * 2);
  viennacl::vector<NumericT> tmp2(vec.size() * 2);

  viennacl::linalg::real_to_complex(vec, tmp, vec.size());
  viennacl::linalg::convolve(circ, tmp, tmp2);
  viennacl::linalg::complex_to_real(tmp2, result, vec.size());

}

typename result_of::size_type< viennacl::vector_base<T> >::type
stride(viennacl::vector_base<T> const & s) { return s.stride(); }

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);

/**
  *  In the following we iterate over all \f$ N \times N \f$ points and apply the five-point stencil directly.
  *  This is done in a straightforward manner for illustration purposes.
  *  Multi-threaded execution via OpenMP can be obtained by uncommenting the pragma below.
  *
  *  Feel free to apply additional optimizations with respect to data access patterns and the like.
  **/

  // 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;
    }
}

 static vcl_size_t size2(viennacl::vector_base<ScalarType> const & /*lhs*/,
                         viennacl::vector_base<ScalarType> const & rhs) { return rhs.size(); }

 static vcl_size_t size1(viennacl::vector_base<ScalarType> const & lhs,
                         viennacl::vector_base<ScalarType> const & /*rhs*/) { return lhs.size(); }