1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
50 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
51 void amg_coarse(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
69 template <
typename InternalType1,
typename InternalType2>
70 void amg_influence(
unsigned int level, InternalType1
const & A, InternalType2 & Pointvector,
amg_tag & tag)
72 typedef typename InternalType1::value_type SparseMatrixType;
73 typedef typename InternalType2::value_type PointVectorType;
74 typedef typename SparseMatrixType::value_type ScalarType;
75 typedef typename SparseMatrixType::value_type ScalarType;
76 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
77 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
84 #pragma omp parallel for private (max,diag_sign) shared (A,Pointvector)
86 for (
unsigned int i=0; i<A[level].size1(); ++i)
89 if (A[level](i,i) < 0)
92 ConstRowIterator row_iter = A[level].begin1();
96 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
98 if (i == (
unsigned int) col_iter.index2())
continue;
100 if (max > *col_iter) max = *col_iter;
102 if (max < *col_iter) max = *col_iter;
110 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
112 unsigned int j = col_iter.index2();
113 if (i == j)
continue;
114 if (diag_sign * (-*col_iter) >= tag.
get_threshold() * (diag_sign * (-max)))
117 Pointvector[level][i]->add_influencing_point(Pointvector[level][j]);
123 std::cout <<
"Influence Matrix: " << std::endl;
124 boost::numeric::ublas::matrix<bool> mat;
125 Pointvector[level].get_influence_matrix(mat);
130 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
132 for (
typename amg_point::iterator iter2 = (*iter)->begin_influencing(); iter2 != (*iter)->end_influencing(); ++iter2)
134 (*iter2)->add_influenced_point(*iter);
139 std::cout <<
"Influence Measures: " << std::endl;
140 boost::numeric::ublas::vector<unsigned int> temp;
141 Pointvector[level].get_influence(temp);
143 std::cout <<
"Point Sorting: " << std::endl;
144 Pointvector[level].get_sorting(temp);
155 template <
typename InternalType1,
typename InternalType2>
158 typedef typename InternalType1::value_type SparseMatrixType;
159 typedef typename InternalType2::value_type PointVectorType;
160 typedef typename SparseMatrixType::value_type ScalarType;
162 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
163 typedef typename SparseMatrixType::iterator2 InternalColIterator;
173 #pragma omp parallel for private (i) shared (Pointvector)
175 for (i=0; i<Pointvector[level].size(); ++i)
176 Pointvector[level][i]->calc_influence();
179 Pointvector[level].sort();
182 while ((c_point = Pointvector[level].get_nextpoint()) != NULL)
185 Pointvector[level].make_cpoint(c_point);
202 Pointvector[level].add_influence(point2,1);
224 #if defined (DEBUG)// or defined (DEBUGBENCH)
225 unsigned int c_points = Pointvector[level].get_cpoints();
226 unsigned int f_points = Pointvector[level].get_fpoints();
227 std::cout <<
"1st pass: Level " << level <<
": ";
228 std::cout <<
"No of C points = " << c_points <<
", ";
229 std::cout <<
"No of F points = " << f_points << std::endl;
233 std::cout <<
"Coarse Points:" << std::endl;
234 boost::numeric::ublas::vector<bool> C;
235 Pointvector[level].get_C(C);
237 std::cout <<
"Fine Points:" << std::endl;
238 boost::numeric::ublas::vector<bool> F;
239 Pointvector[level].get_F(F);
250 template <
typename InternalType1,
typename InternalType2>
253 typedef typename InternalType1::value_type SparseMatrixType;
254 typedef typename InternalType2::value_type PointVectorType;
255 typedef typename SparseMatrixType::value_type ScalarType;
257 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
258 typedef typename SparseMatrixType::iterator2 InternalColIterator;
267 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
293 if ((*iter2)->get_index() == (*iter3)->get_index())
299 else if ((*iter2)->get_index() < (*iter3)->get_index())
339 Pointvector[level].switch_ftoc(point2);
346 std::cout <<
"After 2nd pass:" << std::endl;
347 std::cout <<
"Coarse Points:" << std::endl;
348 boost::numeric::ublas::vector<bool> C;
349 Pointvector[level].get_C(C);
351 std::cout <<
"Fine Points:" << std::endl;
352 boost::numeric::ublas::vector<bool> F;
353 Pointvector[level].get_F(F);
362 std::cout <<
"No C and no F point: ";
363 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
364 if ((*iter)->is_undecided())
365 std::cout << (*iter)->get_index() <<
" ";
366 std::cout << std::endl;
378 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
379 void amg_coarse_rs0(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
381 typedef typename InternalType1::value_type SparseMatrixType;
382 typedef typename InternalType2::value_type PointVectorType;
383 typedef typename SparseMatrixType::value_type ScalarType;
385 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
386 typedef typename SparseMatrixType::iterator2 InternalColIterator;
388 unsigned int total_points;
391 Slicing.slice(level, A, Pointvector);
396 #pragma omp parallel for shared (total_points,Slicing,level)
398 for (
unsigned int i=0; i<Slicing._threads; ++i)
404 Slicing.Offset[i+1][level+1] = Slicing.Pointvector_slice[i][level].get_cpoints();
408 total_points += Slicing.Pointvector_slice[i][level].get_cpoints();
412 if (total_points != 0)
415 #pragma omp parallel for shared (Slicing)
417 for (
unsigned int i=0; i<Slicing._threads; ++i)
420 if (Slicing.Offset[i+1][level+1] == 0)
423 for (
unsigned int j=0; j<Slicing.A_slice[i][level].size1(); ++j)
424 Slicing.Pointvector_slice[i][level].make_cpoint(Slicing.Pointvector_slice[i][level][j]);
425 Slicing.Offset[i+1][level+1] = Slicing.A_slice[i][level].size1();
430 for (
unsigned int i=2; i<=Slicing._threads; ++i)
431 Slicing.Offset[i][level+1] += Slicing.Offset[i-1][level+1];
434 Slicing.join(level, Pointvector);
440 #if defined(DEBUG)// or defined (DEBUGBENCH)
441 for (
unsigned int i=0; i<Slicing._threads; ++i)
443 unsigned int c_points = Slicing.Pointvector_slice[i][level].get_cpoints();
444 unsigned int f_points = Slicing.Pointvector_slice[i][level].get_fpoints();
445 std::cout <<
"Thread " << i <<
": ";
446 std::cout <<
"No of C points = " << c_points <<
", ";
447 std::cout <<
"No of F points = " << f_points << std::endl;
459 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
460 void amg_coarse_rs3(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
462 typedef typename InternalType1::value_type SparseMatrixType;
463 typedef typename InternalType2::value_type PointVectorType;
464 typedef typename SparseMatrixType::value_type ScalarType;
466 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
467 typedef typename SparseMatrixType::iterator2 InternalColIterator;
477 boost::numeric::ublas::vector<unsigned int> Offset = boost::numeric::ublas::vector<unsigned int> (Slicing.Offset.size());
478 for (i=0; i<Slicing.Offset.size(); ++i)
479 Offset[i] = Slicing.Offset[i][level];
482 for (i=0; i<Slicing._threads; ++i)
485 for (j=Offset[i]; j<Offset[i+1]; ++j)
487 point1 = Pointvector[level][j];
511 if ((*iter2)->get_index() == (*iter3)->get_index())
517 else if ((*iter2)->get_index() < (*iter3)->get_index())
564 Pointvector[level].switch_ftoc(point2);
566 for (
unsigned int j=i+1; j<=Slicing._threads; ++j)
567 Slicing.Offset[j][level+1]++;
576 std::cout <<
"After 3rd pass:" << std::endl;
577 std::cout <<
"Coarse Points:" << std::endl;
578 boost::numeric::ublas::vector<bool> C;
579 Pointvector[level].get_C(C);
581 std::cout <<
"Fine Points:" << std::endl;
582 boost::numeric::ublas::vector<bool> F;
583 Pointvector[level].get_F(F);
593 std::cout <<
"No C and no F point: ";
594 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
595 if ((*iter)->is_undecided())
596 std::cout << i <<
" ";
597 std::cout << std::endl;
609 template <
typename InternalType1,
typename InternalType2>
612 typedef typename InternalType1::value_type SparseMatrixType;
613 typedef typename InternalType2::value_type PointVectorType;
614 typedef typename SparseMatrixType::value_type ScalarType;
616 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
617 typedef typename SparseMatrixType::iterator2 InternalColIterator;
624 if (A[level].
size1() == 1)
return;
629 #pragma omp parallel for private (x,y,diag) shared (A)
631 for (x=0; x<A[level].size1(); ++x)
633 InternalRowIterator row_iter = A[level].begin1();
635 diag = A[level](x,x);
636 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
638 y = col_iter.index2();
639 if (y == x || (std::abs(*col_iter) >= tag.
get_threshold()*pow(0.5,level-1) * sqrt(std::abs(diag*A[level](y,y)))))
648 std::cout <<
"Neighborhoods:" << std::endl;
649 boost::numeric::ublas::matrix<bool> mat;
650 Pointvector[level].get_influence_matrix(mat);
655 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
679 std::cout <<
"After aggregation:" << std::endl;
680 std::cout <<
"Coarse Points:" << std::endl;
681 boost::numeric::ublas::vector<bool> C;
682 Pointvector[level].get_C(C);
684 std::cout <<
"Fine Points:" << std::endl;
685 boost::numeric::ublas::vector<bool> F;
686 Pointvector[level].get_F(F);
688 std::cout <<
"Aggregates:" << std::endl;