1 #ifndef VIENNACL_LINALG_DETAIL_SPAI_SPAI_HPP
2 #define VIENNACL_LINALG_DETAIL_SPAI_SPAI_HPP
43 #include "boost/numeric/ublas/vector.hpp"
44 #include "boost/numeric/ublas/matrix.hpp"
45 #include "boost/numeric/ublas/matrix_proxy.hpp"
46 #include "boost/numeric/ublas/vector_proxy.hpp"
47 #include "boost/numeric/ublas/storage.hpp"
48 #include "boost/numeric/ublas/io.hpp"
49 #include "boost/numeric/ublas/lu.hpp"
50 #include "boost/numeric/ublas/triangular.hpp"
51 #include "boost/numeric/ublas/matrix_expression.hpp"
68 #define VIENNACL_SPAI_K_b 20
80 template<
typename SparseVectorType>
82 for(
typename SparseVectorType::const_iterator vec_it = v.begin(); vec_it!= v.end(); ++vec_it){
83 std::cout<<
"[ "<<vec_it->first<<
" ]:"<<vec_it->second<<std::endl;
86 template<
typename DenseMatrixType>
88 for(
int i = 0; i < m.size2(); ++i){
89 for(
int j = 0; j < m.size1(); ++j){
90 std::cout<<m(j, i)<<
" ";
101 template<
typename SparseVectorType,
typename ScalarType>
103 for(
typename SparseVectorType::const_iterator v_it = v.begin(); v_it != v.end(); ++v_it){
104 res_v[v_it->first] += b*v_it->second;
114 template<
typename SparseVectorType,
typename ScalarType>
116 const unsigned int ind, SparseVectorType& res){
117 for(
typename SparseVectorType::const_iterator v_it = v.begin(); v_it != v.end(); ++v_it){
120 res[ind] -=
static_cast<ScalarType
>(1);
129 template<
typename SparseVectorType>
130 void build_index_set(
const std::vector<SparseVectorType>& A_v_c,
const SparseVectorType& v, std::vector<unsigned int>& J,
131 std::vector<unsigned int>& I){
146 template<
typename SparseMatrixType,
typename DenseMatrixType,
typename SparseVectorType,
typename VectorType>
148 const std::vector<SparseVectorType>& A_v_c,
149 const std::vector<SparseVectorType>& M_v,
150 std::vector<std::vector<unsigned int> >& g_I,
151 std::vector<std::vector<unsigned int> >& g_J,
152 std::vector<DenseMatrixType>& g_A_I_J,
153 std::vector<VectorType>& g_b_v){
155 #pragma omp parallel for
157 for(std::size_t i = 0; i < M_v.size(); ++i){
172 template<
typename SparseVectorType>
174 const std::vector<SparseVectorType>& M_v,
175 std::vector<std::vector<unsigned int> >& g_J,
176 std::vector<std::vector<unsigned int> >& g_I)
179 #pragma omp parallel for
181 for(std::size_t i = 0; i < M_v.size(); ++i){
198 template<
typename ScalarType,
unsigned int MAT_ALIGNMENT,
typename SparseVectorType>
200 const std::vector<SparseVectorType>& A_v_c,
201 const std::vector<SparseVectorType>& M_v,
202 std::vector<cl_uint> g_is_update,
203 std::vector<std::vector<unsigned int> >& g_I,
204 std::vector<std::vector<unsigned int> >& g_J,
207 const unsigned int cur_iter)
212 block_assembly(A, g_J, g_I, g_A_I_J, g_is_update, is_empty_block, cur_iter);
213 block_qr<ScalarType>(g_I, g_J, g_A_I_J, g_bv, g_is_update, cur_iter);
227 template<
typename ScalarType,
typename SparseVectorType>
229 unsigned int start_m_ind,
230 const std::vector<unsigned int> & J,
231 SparseVectorType & m)
233 unsigned int cnt = 0;
234 for (std::size_t i = 0; i < J.size(); ++i) {
235 m[J[i]] = m_in[start_m_ind + cnt++];
254 template<
typename SparseVectorType,
typename ScalarType>
256 std::vector<SparseVectorType> & M_v,
257 std::vector<std::vector<unsigned int> >& g_I,
258 std::vector<std::vector<unsigned int> > & g_J,
261 std::vector<SparseVectorType> & g_res,
262 std::vector<cl_uint> & g_is_update,
264 const unsigned int cur_iter){
265 unsigned int y_sz, m_sz;
266 std::vector<cl_uint> y_inds(M_v.size() + 1,
static_cast<cl_uint
>(0));
267 std::vector<cl_uint> m_inds(M_v.size() + 1,
static_cast<cl_uint
>(0));
272 std::vector<ScalarType> y_v(y_sz, static_cast<ScalarType>(0));
273 for(std::size_t i = 0; i < M_v.size(); ++i){
274 for(std::size_t j = 0; j < g_I[i].size(); ++j){
276 y_v[y_inds[i] + j] =
static_cast<ScalarType
>(1.0);
281 std::vector<ScalarType> m_v(m_sz, static_cast<cl_uint>(0));
297 static_cast<unsigned int>(
sizeof(ScalarType)*y_v.size()),
300 static_cast<unsigned int>(
sizeof(ScalarType)*m_v.size()),
303 static_cast<unsigned int>(
sizeof(cl_uint)*(g_I.size() + 1)),
306 static_cast<unsigned int>(
sizeof(cl_uint)*(g_is_update.size())),
313 g_A_I_J_vcl.
handle1(), g_is_update_vcl,
315 static_cast<unsigned int>(M_v.size())));
318 m_v_vcl.
handle(), CL_TRUE, 0,
319 sizeof(ScalarType)*(m_v.size()),
320 &(m_v[0]), 0, NULL, NULL);
324 for(std::size_t i = 0; i < M_v.size(); ++i){
329 compute_spai_residual<SparseVectorType, ScalarType>(A_v_c, M_v[i],
static_cast<unsigned int>(i), g_res[i]);
330 ScalarType res_norm = 0;
354 template<
typename SparseVectorType,
typename DenseMatrixType,
typename VectorType>
356 std::vector<DenseMatrixType>& g_R,
357 std::vector<VectorType>& g_b_v,
358 std::vector<std::vector<unsigned int> >& g_I,
359 std::vector<std::vector<unsigned int> >& g_J,
360 std::vector<SparseVectorType>& g_res,
361 std::vector<bool>& g_is_update,
362 std::vector<SparseVectorType>& M_v,
364 typedef typename DenseMatrixType::value_type ScalarType;
367 #pragma omp parallel for
369 for(std::size_t i = 0; i < M_v.size(); ++i){
371 VectorType y = boost::numeric::ublas::zero_vector<ScalarType>(g_I[i].size());
373 projectI<VectorType, ScalarType>(g_I[i], y,
static_cast<unsigned int>(tag.
getBegInd() + i));
375 VectorType m_new = boost::numeric::ublas::zero_vector<ScalarType>(g_R[i].size2());
379 compute_spai_residual<SparseVectorType, ScalarType>(A_v_c, M_v[i],
static_cast<unsigned int>(tag.
getBegInd() + i), g_res[i]);
380 ScalarType res_norm = 0;
391 template<
typename VectorType>
394 for(
unsigned int i = 0; i < parallel_is_update.size(); ++i){
395 if(parallel_is_update[i])
407 template<
typename SparseMatrixType,
typename SparseVectorType>
409 for(
typename SparseMatrixType::const_iterator1 row_it = M_in.begin1(); row_it!= M_in.end1(); ++row_it){
411 for(
typename SparseMatrixType::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it){
412 M_v[
static_cast<unsigned int>(col_it.index2())][static_cast<unsigned int>(col_it.index1())] = *col_it;
419 template<
typename SparseMatrixType,
typename SparseVectorType>
421 for(
typename SparseMatrixType::const_iterator1 row_it = M_in.begin1(); row_it!= M_in.end1(); ++row_it){
422 for(
typename SparseMatrixType::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it){
423 M_v[
static_cast<unsigned int>(col_it.index1())][static_cast<unsigned int>(col_it.index2())] = *col_it;
432 void write_set_to_array(
const std::vector<std::vector<unsigned int> >& ind_set, std::vector<cl_uint>& a){
433 unsigned int cnt = 0;
435 for(
size_t i = 0; i < ind_set.size(); ++i){
436 for(
size_t j = 0; j < ind_set[i].size(); ++j){
437 a[cnt++] =
static_cast<cl_uint
>(ind_set[i][j]);
454 template<
typename ScalarType,
unsigned int MAT_ALIGNMENT>
456 const std::vector<std::vector<unsigned int> >& g_I,
458 std::vector<cl_uint>& g_is_update,
459 bool& is_empty_block,
460 const unsigned int cur_iter){
462 unsigned int sz_I, sz_J, sz_blocks;
463 std::vector<cl_uint> matrix_dims(g_I.size()*2,
static_cast<cl_uint
>(0));
464 std::vector<cl_uint> i_ind(g_I.size() + 1,
static_cast<cl_uint
>(0));
465 std::vector<cl_uint> j_ind(g_I.size() + 1,
static_cast<cl_uint
>(0));
466 std::vector<cl_uint> blocks_ind(g_I.size() + 1,
static_cast<cl_uint
>(0));
473 std::vector<cl_uint> I_set(sz_I, static_cast<cl_uint>(0));
475 std::vector<cl_uint> J_set(sz_J, static_cast<cl_uint>(0));
481 if(I_set.size() > 0 && J_set.size() > 0){
483 std::vector<ScalarType> con_A_I_J(sz_blocks, static_cast<ScalarType>(0));
498 static_cast<unsigned int>(
sizeof(ScalarType)*(sz_blocks)),
502 static_cast<unsigned int>(
sizeof(cl_uint)*2*static_cast<cl_uint>(g_I.size())),
506 static_cast<unsigned int>(
sizeof(cl_uint)*(g_I.size() + 1)),
510 static_cast<unsigned int>(
sizeof(cl_uint)*sz_I),
514 static_cast<unsigned int>(
sizeof(cl_uint)*sz_J),
518 static_cast<unsigned int>(
sizeof(cl_uint)*(g_I.size() + 1)),
522 static_cast<unsigned int>(
sizeof(cl_uint)*(g_I.size() + 1)),
526 static_cast<unsigned int>(
sizeof(cl_uint)*g_is_update.size()),
536 static_cast<unsigned int>(g_I.size())));
537 is_empty_block =
false;
539 is_empty_block =
true;
550 template<
typename SparseMatrixType,
typename SparseVectorType>
556 for(
unsigned int i = 0; i < M_v.size(); ++i){
557 for(
typename SparseVectorType::const_iterator vec_it = M_v[i].begin(); vec_it!=M_v[i].end(); ++vec_it){
558 M(vec_it->first, i) = vec_it->second;
564 for(
unsigned int i = 0; i < M_v.size(); ++i){
565 for(
typename SparseVectorType::const_iterator vec_it = M_v[i].begin(); vec_it!=M_v[i].end(); ++vec_it){
566 M(i, vec_it->first) = vec_it->second;
576 template<
typename MatrixType>
578 typedef typename MatrixType::value_type ScalarType;
579 std::vector<std::map<size_t, ScalarType> > temp_A(A_in.size2());
580 A.resize(A_in.size2(), A_in.size1(),
false);
582 for (
typename MatrixType::const_iterator1 row_it = A_in.begin1();
583 row_it != A_in.end1();
586 for (
typename MatrixType::const_iterator2 col_it = row_it.begin();
587 col_it != row_it.end();
590 temp_A[col_it.index2()][col_it.index1()] = *col_it;
594 for (
size_t i=0; i<temp_A.size(); ++i)
596 for (
typename std::map<size_t, ScalarType>::const_iterator it = temp_A[i].begin();
597 it != temp_A[i].end();
599 A(i, it->first) = it->second;
619 template <
typename MatrixType>
621 typedef typename MatrixType::value_type ScalarType;
622 typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
624 typedef typename boost::numeric::ublas::matrix<ScalarType> DenseMatrixType;
626 unsigned int cur_iter = 0;
629 std::vector<SparseVectorType> A_v_c(M.size2());
630 std::vector<SparseVectorType> M_v(M.size2());
636 go_on = (tag.
getEndInd() <
static_cast<long>(M.size2()));
640 std::vector<bool> g_is_update(l_sz,
true);
647 std::vector<SparseVectorType> l_M_v(l_sz);
658 std::vector<DenseMatrixType> g_A_I_J(l_sz);
660 std::vector<VectorType> g_b_v(l_sz);
662 std::vector<SparseVectorType> g_res(l_sz);
664 std::vector<std::vector<unsigned int> > g_I(l_sz);
666 std::vector<std::vector<unsigned int> > g_J(l_sz);
670 if(cur_iter == 0)
block_set_up(A, A_v_c, l_M_v, g_I, g_J, g_A_I_J, g_b_v);
671 else block_update(A, A_v_c, g_res, g_is_update, g_I, g_J, g_b_v, g_A_I_J, tag);
687 M.resize(M.size1(), M.size2(),
false);
700 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
702 const boost::numeric::ublas::compressed_matrix<ScalarType>& cpu_A,
703 boost::numeric::ublas::compressed_matrix<ScalarType>& cpu_M,
706 typedef typename boost::numeric::ublas::vector<ScalarType> VectorType;
708 typedef typename boost::numeric::ublas::matrix<ScalarType> DenseMatrixType;
709 typedef typename boost::numeric::ublas::compressed_matrix<ScalarType> CPUMatrixType;
712 unsigned int cur_iter = 0;
713 std::vector<cl_uint> g_is_update(cpu_M.size2(),
static_cast<cl_uint
>(1));
716 std::vector<SparseVectorType> A_v_c(cpu_M.size2());
717 std::vector<SparseVectorType> M_v(cpu_M.size2());
720 std::vector<SparseVectorType> g_res(cpu_M.size2());
721 std::vector<std::vector<unsigned int> > g_I(cpu_M.size2());
722 std::vector<std::vector<unsigned int> > g_J(cpu_M.size2());
732 if(cur_iter == 0)
block_set_up(A, A_v_c, M_v, g_is_update, g_I, g_J, g_A_I_J_vcl, g_bv_vcl, cur_iter);
733 else block_update(A, A_v_c, g_is_update, g_res, g_J, g_I, g_A_I_J_vcl, g_bv_vcl, tag, cur_iter);
738 least_square_solve<SparseVectorType, ScalarType>(A_v_c, M_v, g_I, g_J, g_A_I_J_vcl, g_bv_vcl, g_res, g_is_update, tag, cur_iter);
740 if(tag.getIsStatic())
break;
743 cpu_M.resize(cpu_M.size1(), cpu_M.size2(),
false);
746 M.
resize(static_cast<unsigned int>(cpu_M.size1()), static_cast<unsigned int>(cpu_M.size2()));