1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_QR_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_QR_HPP
36 #include "boost/numeric/ublas/vector.hpp"
37 #include "boost/numeric/ublas/matrix.hpp"
38 #include "boost/numeric/ublas/matrix_proxy.hpp"
39 #include "boost/numeric/ublas/storage.hpp"
40 #include "boost/numeric/ublas/io.hpp"
41 #include "boost/numeric/ublas/matrix_expression.hpp"
42 #include "boost/numeric/ublas/detail/matrix_assign.hpp"
65 template<
typename T,
typename InputIterator>
66 void Print(std::ostream& ostr, InputIterator it_begin, InputIterator it_end){
68 std::string delimiters =
" ";
69 std::copy(it_begin, it_end, std::ostream_iterator<T>(ostr, delimiters.c_str()));
73 template<
typename VectorType,
typename MatrixType>
74 void write_to_block(VectorType& con_A_I_J,
unsigned int start_ind,
const std::vector<unsigned int>& I,
const std::vector<unsigned int>& J, MatrixType& m){
75 m.resize(I.size(), J.size(),
false);
76 for(
size_t i = 0; i < J.size(); ++i){
77 for(
size_t j = 0; j < I.size(); ++j){
78 m(j,i) = con_A_I_J[start_ind + i*I.size() + j];
83 template<
typename VectorType>
85 const std::vector<std::vector<unsigned int> >& g_I,
const std::vector<std::vector<unsigned int> >& g_J){
86 typedef typename VectorType::value_type ScalarType;
87 std::vector<boost::numeric::ublas::matrix<ScalarType> > com_A_I_J(g_I.size());
88 for(
size_t i = 0; i < g_I.size(); ++i){
89 write_to_block( con_A_I_J, blocks_ind[i], g_I[i], g_J[i], com_A_I_J[i]);
90 std::cout<<com_A_I_J[i]<<std::endl;
93 template<
typename VectorType>
94 void print_continious_vector(VectorType& con_v, std::vector<cl_uint>& block_ind,
const std::vector<std::vector<unsigned int> >& g_J){
95 typedef typename VectorType::value_type ScalarType;
96 std::vector<boost::numeric::ublas::vector<ScalarType> > com_v(g_J.size());
98 for(
size_t i = 0; i < g_J.size(); ++i){
99 com_v[i].resize(g_J[i].
size());
100 for(
size_t j = 0; j < g_J[i].size(); ++j){
101 com_v[i](j) = con_v[block_ind[i] + j];
103 std::cout<<com_v[i]<<std::endl;
115 void compute_blocks_size(
const std::vector<std::vector<unsigned int> >& g_I,
const std::vector<std::vector<unsigned int> >& g_J,
116 unsigned int& sz, std::vector<cl_uint>& blocks_ind, std::vector<cl_uint>& matrix_dims){
118 for(
size_t i = 0; i < g_I.size(); ++i){
119 sz +=
static_cast<unsigned int>(g_I[i].size()*g_J[i].size());
120 matrix_dims[2*i] =
static_cast<cl_uint
>(g_I[i].size());
121 matrix_dims[2*i + 1] =
static_cast<cl_uint
>(g_J[i].size());
122 blocks_ind[i+1] = blocks_ind[i] +
static_cast<cl_uint
>(g_I[i].size()*g_J[i].size());
130 void get_size(
const std::vector<std::vector<unsigned int> >& inds,
unsigned int&
size){
132 for (
size_t i = 0; i < inds.size(); ++i) {
133 size +=
static_cast<unsigned int>(inds[i].size());
141 void init_start_inds(
const std::vector<std::vector<unsigned int> >& inds, std::vector<cl_uint>& start_inds){
142 for(
size_t i = 0; i < inds.size(); ++i){
143 start_inds[i+1] = start_inds[i] +
static_cast<cl_uint
>(inds[i].size());
153 template<
typename MatrixType,
typename ScalarType>
154 void dot_prod(
const MatrixType& A,
unsigned int beg_ind, ScalarType& res){
155 res =
static_cast<ScalarType
>(0);
156 for(
size_t i = beg_ind; i < A.size1(); ++i){
157 res += A(i, beg_ind-1)*A(i, beg_ind-1);
167 template<
typename MatrixType,
typename VectorType,
typename ScalarType>
168 void custom_inner_prod(
const MatrixType& A,
const VectorType& v,
unsigned int col_ind,
unsigned int start_ind, ScalarType& res){
169 res =
static_cast<ScalarType
>(0);
170 for(
unsigned int i = start_ind; i < static_cast<unsigned int>(A.size1()); ++i){
171 res += A(i, col_ind)*v(i);
180 template<
typename MatrixType,
typename VectorType>
181 void copy_vector(
const MatrixType & A, VectorType & v,
const unsigned int beg_ind){
182 for(
unsigned int i = beg_ind; i < static_cast<unsigned int>(A.size1()); ++i){
183 v(i) = A( i, beg_ind-1);
194 template<
typename MatrixType,
typename VectorType,
typename ScalarType>
201 v(j) =
static_cast<ScalarType
>(1.0);
206 mu = std::sqrt(A(j,j)*A(j, j) + sg);
210 v(j) = -sg/(A(j, j) + mu);
212 b = 2*(v(j)*v(j))/(sg + v(j)*v(j));
222 template<
typename MatrixType,
typename VectorType,
typename ScalarType>
225 ScalarType in_prod_res;
226 for(
unsigned int i = iter_cnt; i < static_cast<unsigned int>(A.size2()); ++i){
229 for(
unsigned int j = iter_cnt; j < static_cast<unsigned int>(A.size1()); ++j){
230 A(j, i) -= b*in_prod_res*v(j);
240 template<
typename MatrixType,
typename VectorType>
242 for(
unsigned int i = ind; i < static_cast<unsigned int>(A.size1()); ++i){
253 template<
typename MatrixType,
typename VectorType>
255 typedef typename MatrixType::value_type ScalarType;
256 if((R.size1() > 0) && (R.size2() > 0)){
257 VectorType v = (VectorType)boost::numeric::ublas::zero_vector<ScalarType>(R.size1());
258 b_v = (VectorType)boost::numeric::ublas::zero_vector<ScalarType>(R.size2());
259 for(
unsigned int i = 0; i < static_cast<unsigned int>(R.size2()); ++i){
273 std::ifstream ifs(file_name.c_str(), std::ifstream::in);
276 std::cerr <<
"WARNING: Cannot open file " << file_name << std::endl;
279 std::ostringstream ost;
280 while (std::getline(ifs, line)) {
281 ost<<line<<std::endl;
283 kernel_source = ost.str();
290 void get_max_block_size(
const std::vector<std::vector<unsigned int> >& inds,
unsigned int& max_size){
292 for(
unsigned int i = 0; i < inds.size(); ++i){
293 if(inds[i].
size() > max_size){
294 max_size =
static_cast<unsigned int>(inds[i].size());
305 template<
typename MatrixType,
typename VectorType,
typename ScalarType>
306 void custom_dot_prod(
const MatrixType& A,
const VectorType& v,
unsigned int ind, ScalarType& res){
307 res =
static_cast<ScalarType
>(0);
308 for(
unsigned int j = ind; j < A.size1(); ++j){
312 res += A(j, ind)*v(j);
322 template<
typename MatrixType,
typename VectorType>
324 typedef typename MatrixType::value_type ScalarType;
325 ScalarType inn_prod =
static_cast<ScalarType
>(0);
326 for(
size_t i = 0; i < R.size2(); ++i){
328 for(
size_t j = i; j < R.size1(); ++j){
330 y(j) -= b_v(i)*inn_prod;
333 y(j) -= b_v(i)*inn_prod*R(j,i);
344 template<
typename MatrixType,
typename VectorType>
347 for(
size_t i = 0; i < A.size2(); ++i){
348 tmp_v = (VectorType)column(A,i);
364 template<
typename ScalarType>
365 void block_qr(std::vector<std::vector<unsigned int> >& g_I,
366 std::vector<std::vector<unsigned int> >& g_J,
369 std::vector<cl_uint>& g_is_update,
370 const unsigned int cur_iter){
372 unsigned int bv_size;
376 unsigned int local_r_n, local_c_n;
384 std::vector<cl_uint> start_bv_inds(g_I.size() + 1, 0);
385 std::vector<cl_uint> start_v_inds(g_I.size() + 1, 0);
389 std::vector<ScalarType> b_v(bv_size, static_cast<ScalarType>(0));
390 std::vector<ScalarType> v(v_size, static_cast<ScalarType>(0));
405 static_cast<unsigned int>(
sizeof(ScalarType)*bv_size),
409 static_cast<unsigned int>(
sizeof(ScalarType)*v_size),
413 static_cast<unsigned int>(
sizeof(cl_uint)*g_I.size()),
414 &(start_bv_inds[0]));
417 static_cast<unsigned int>(
sizeof(cl_uint)*g_I.size()),
420 static_cast<unsigned int>(
sizeof(cl_uint)*g_is_update.size()),
431 static_cast<cl_uint
>(g_I.size())));