C++ matrix::internal_size2方法代码示例

double luT(
	viennacl::matrix<T> &vclX,
	viennacl::vector_base<T> &vclD) {

	double logdet=-99.9;

	viennacl::linalg::lu_factorize(vclX);

	// pointer to the actual diagonal
	viennacl::vector_base<T> diagOfVar(
		vclX.handle(), vclX.size1(), 0, 
		vclX.internal_size2() + 1);

		// compute log determinant
	vclD = viennacl::linalg::element_log(diagOfVar);
	logdet = viennacl::linalg::sum(vclD);
// OPERATION_UNARY_LOG_TYPE 	
		//http://viennacl.sourceforge.net/doc/scheduler_8cpp-example.html#a11

		// put the diagonals in D, and 1's on the diagonal of L
	vclD = diagOfVar;

		//diagOfVar = T(1); // problem here
	
	return(logdet);
}

void init_random(viennacl::matrix<T, F> & M)
{
    std::vector<T> cM(M.internal_size());
    for (std::size_t i = 0; i < M.size1(); ++i)
        for (std::size_t j = 0; j < M.size2(); ++j)
            cM[F::mem_index(i, j, M.internal_size1(), M.internal_size2())] = T(rand())/T(RAND_MAX);
    viennacl::fast_copy(&cM[0],&cM[0] + cM.size(),M);
}

      bool householder_c(viennacl::matrix<SCALARTYPE, row_major, ALIGNMENT>& A,
                          viennacl::matrix<SCALARTYPE, row_major, ALIGNMENT>& Q,
                          viennacl::vector<SCALARTYPE, ALIGNMENT>& D,
                          vcl_size_t row_start, vcl_size_t col_start)
      {
        viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(A).context());

        if (row_start + 1 >= A.size1())
          return false;

        prepare_householder_vector(A, D, A.size1(), row_start, col_start, row_start, true);

        {
          viennacl::ocl::kernel& kernel = ctx.get_kernel(viennacl::linalg::opencl::kernels::svd<SCALARTYPE>::program_name(), SVD_HOUSEHOLDER_UPDATE_A_LEFT_KERNEL);

          viennacl::ocl::enqueue(kernel(
                                        A,
                                        D,
                                        static_cast<cl_uint>(row_start),
                                        static_cast<cl_uint>(col_start),
                                        static_cast<cl_uint>(A.size1()),
                                        static_cast<cl_uint>(A.size2()),
                                        static_cast<cl_uint>(A.internal_size2()),
                                        viennacl::ocl::local_mem(static_cast<cl_uint>(128 * sizeof(SCALARTYPE)))
                                ));

        }

        {
          viennacl::ocl::kernel& kernel = ctx.get_kernel(viennacl::linalg::opencl::kernels::svd<SCALARTYPE>::program_name(), SVD_HOUSEHOLDER_UPDATE_QL_KERNEL);

          viennacl::ocl::enqueue(kernel(
                                        Q,
                                        D,
                                        static_cast<cl_uint>(A.size1()),
                                      //  static_cast<cl_uint>(A.size2()),
                                        static_cast<cl_uint>(Q.internal_size2()),
                                        viennacl::ocl::local_mem(static_cast<cl_uint>(128 * sizeof(SCALARTYPE)))
                                ));

        }

        return true;
      }

NumericT diff(std::vector<std::vector<NumericT> > const & A1, viennacl::matrix<NumericT> const & A2)
{
  std::vector<NumericT> host_values(A2.internal_size());
  for (std::size_t i=0; i<A2.size1(); ++i)
    for (std::size_t j=0; j<A2.size2(); ++j)
      host_values[i*A2.internal_size2() + j] = A1[i][j];

  std::vector<NumericT> device_values(A2.internal_size());
  viennacl::fast_copy(A2, &device_values[0]);
  viennacl::vector<NumericT> vcl_device_values(A2.internal_size());  // workaround to avoid code duplication
  viennacl::copy(device_values, vcl_device_values);

  return diff(host_values, vcl_device_values);
}

    void nmf(viennacl::matrix<ScalarType> const & v,
             viennacl::matrix<ScalarType> & w,
             viennacl::matrix<ScalarType> & h,
             std::size_t k,
             ScalarType eps = 0.000001,
             std::size_t max_iter = 10000,
             std::size_t check_diff_every_step = 100)
    {
      viennacl::linalg::kernels::nmf<ScalarType, 1>::init();
      
      w.resize(v.size1(), k);
      h.resize(k, v.size2());

      std::vector<ScalarType> stl_w(w.internal_size1() * w.internal_size2());
      std::vector<ScalarType> stl_h(h.internal_size1() * h.internal_size2());

      for (std::size_t j = 0; j < stl_w.size(); j++)
          stl_w[j] = static_cast<ScalarType>(rand()) / RAND_MAX;

      for (std::size_t j = 0; j < stl_h.size(); j++)
          stl_h[j] = static_cast<ScalarType>(rand()) / RAND_MAX;

      viennacl::matrix<ScalarType> wn(v.size1(), k);
      viennacl::matrix<ScalarType> wd(v.size1(), k);
      viennacl::matrix<ScalarType> wtmp(v.size1(), v.size2());

      viennacl::matrix<ScalarType> hn(k, v.size2());
      viennacl::matrix<ScalarType> hd(k, v.size2());
      viennacl::matrix<ScalarType> htmp(k, k);

      viennacl::matrix<ScalarType> appr(v.size1(), v.size2());
      viennacl::vector<ScalarType> diff(v.size1() * v.size2());

      viennacl::fast_copy(&stl_w[0], &stl_w[0] + stl_w.size(), w);
      viennacl::fast_copy(&stl_h[0], &stl_h[0] + stl_h.size(), h);

      ScalarType last_diff = 0.0f;


      
      for (std::size_t i = 0; i < max_iter; i++)
      {
        {
          hn = viennacl::linalg::prod(trans(w), v);
          htmp = viennacl::linalg::prod(trans(w), w);
          hd = viennacl::linalg::prod(htmp, h);

          viennacl::ocl::kernel & mul_div_kernel = viennacl::ocl::get_kernel(viennacl::linalg::kernels::nmf<ScalarType, 1>::program_name(), 
                                                                             NMF_MUL_DIV_KERNEL);
          viennacl::ocl::enqueue(mul_div_kernel(h, hn, hd, cl_uint(stl_h.size())));
        }
        {
          wn = viennacl::linalg::prod(v, trans(h));
          wtmp = viennacl::linalg::prod(w, h);
          wd = viennacl::linalg::prod(wtmp, trans(h));

          viennacl::ocl::kernel & mul_div_kernel = viennacl::ocl::get_kernel(viennacl::linalg::kernels::nmf<ScalarType, 1>::program_name(), 
                                                                             NMF_MUL_DIV_KERNEL);
          
          viennacl::ocl::enqueue(mul_div_kernel(w, wn, wd, cl_uint(stl_w.size())));
        }

        if (i % check_diff_every_step == 0)
        {
          appr = viennacl::linalg::prod(w, h);

         viennacl::ocl::kernel & sub_kernel = viennacl::ocl::get_kernel(viennacl::linalg::kernels::nmf<ScalarType, 1>::program_name(), 
                                                                        NMF_SUB_KERNEL);
          //this is a cheat. i.e save difference of two matrix into vector to get norm_2
          viennacl::ocl::enqueue(sub_kernel(appr, v, diff, cl_uint(v.size1() * v.size2())));
          ScalarType diff_val = viennacl::linalg::norm_2(diff);

          if((diff_val < eps) || (fabs(diff_val - last_diff) < eps))
          {
              //std::cout << "Breaked at diff - " << diff_val << "\n";
              break;
          }

          last_diff = diff_val;

          //printf("Iteration #%lu - %.5f \n", i, diff_val);
        }
      }
      
      
    }

 viennacl::vector_range<viennacl::vector_base<T> > sharedCol(){
     // viennacl::vector_base<T> tmp(ptr_matrix->handle(), ptr_matrix->internal_size(), 0, 1);
     
     // std::cout << "returning column" << std::endl;
     viennacl::vector_base<T> tmp(ptr_matrix->handle(), ptr_matrix->size1(), begin, ptr_matrix->internal_size2());
     // std::cout << "got column" << std::endl;
     
     viennacl::vector_range<viennacl::vector_base<T> > v_sub(tmp, r);
     // std::cout << "got range" << std::endl;
     
     return v_sub;
 }