1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
29 #ifdef VIENNACL_WITH_OPENMP
51 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
52 void amg_coarse(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
70 template <
typename InternalType1,
typename InternalType2>
71 void amg_influence(
unsigned int level, InternalType1
const & A, InternalType2 & Pointvector,
amg_tag & tag)
73 typedef typename InternalType1::value_type SparseMatrixType;
74 typedef typename InternalType2::value_type PointVectorType;
75 typedef typename SparseMatrixType::value_type ScalarType;
76 typedef typename SparseMatrixType::value_type ScalarType;
77 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
78 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
84 #ifdef VIENNACL_WITH_OPENMP
85 #pragma omp parallel for private (max,diag_sign)
87 for (
long i=0; i<static_cast<long>(A[level].size1()); ++i)
90 if (A[level](i,i) < 0)
93 ConstRowIterator row_iter = A[level].begin1();
97 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
99 if (i == (
unsigned int) col_iter.index2())
continue;
101 if (max > *col_iter) max = *col_iter;
103 if (max < *col_iter) max = *col_iter;
111 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
113 unsigned int j =
static_cast<unsigned int>(col_iter.index2());
114 if (i == j)
continue;
115 if (diag_sign * (-*col_iter) >= tag.
get_threshold() * (diag_sign * (-max)))
118 Pointvector[level][i]->add_influencing_point(Pointvector[level][j]);
123 #ifdef VIENNACL_AMG_DEBUG
124 std::cout <<
"Influence Matrix: " << std::endl;
125 boost::numeric::ublas::matrix<bool> mat;
126 Pointvector[level].get_influence_matrix(mat);
131 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
133 for (
typename amg_point::iterator iter2 = (*iter)->begin_influencing(); iter2 != (*iter)->end_influencing(); ++iter2)
135 (*iter2)->add_influenced_point(*iter);
139 #ifdef VIENNACL_AMG_DEBUG
140 std::cout <<
"Influence Measures: " << std::endl;
141 boost::numeric::ublas::vector<unsigned int> temp;
142 Pointvector[level].get_influence(temp);
144 std::cout <<
"Point Sorting: " << std::endl;
145 Pointvector[level].get_sorting(temp);
156 template <
typename InternalType1,
typename InternalType2>
166 #ifdef VIENNACL_WITH_OPENMP
167 #pragma omp parallel for private (i)
169 for (i=0; i<static_cast<long>(Pointvector[level].size()); ++i)
170 Pointvector[level][i]->calc_influence();
173 Pointvector[level].sort();
176 while ((c_point = Pointvector[level].get_nextpoint()) != NULL)
179 Pointvector[level].make_cpoint(c_point);
196 Pointvector[level].add_influence(point2,1);
218 #if defined (VIENNACL_AMG_DEBUG)// or defined (VIENNACL_AMG_DEBUGBENCH)
219 unsigned int c_points = Pointvector[level].get_cpoints();
220 unsigned int f_points = Pointvector[level].get_fpoints();
221 std::cout <<
"1st pass: Level " << level <<
": ";
222 std::cout <<
"No of C points = " << c_points <<
", ";
223 std::cout <<
"No of F points = " << f_points << std::endl;
226 #ifdef VIENNACL_AMG_DEBUG
227 std::cout <<
"Coarse Points:" << std::endl;
228 boost::numeric::ublas::vector<bool> C;
229 Pointvector[level].get_C(C);
231 std::cout <<
"Fine Points:" << std::endl;
232 boost::numeric::ublas::vector<bool> F;
233 Pointvector[level].get_F(F);
244 template <
typename InternalType1,
typename InternalType2>
247 typedef typename InternalType2::value_type PointVectorType;
256 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
282 if ((*iter2)->get_index() == (*iter3)->get_index())
288 else if ((*iter2)->get_index() < (*iter3)->get_index())
328 Pointvector[level].switch_ftoc(point2);
334 #ifdef VIENNACL_AMG_DEBUG
335 std::cout <<
"After 2nd pass:" << std::endl;
336 std::cout <<
"Coarse Points:" << std::endl;
337 boost::numeric::ublas::vector<bool> C;
338 Pointvector[level].get_C(C);
340 std::cout <<
"Fine Points:" << std::endl;
341 boost::numeric::ublas::vector<bool> F;
342 Pointvector[level].get_F(F);
346 #ifdef VIENNACL_AMG_DEBUG
347 #ifdef VIENNACL_WITH_OPENMP
351 std::cout <<
"No C and no F point: ";
352 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
353 if ((*iter)->is_undecided())
354 std::cout << (*iter)->get_index() <<
" ";
355 std::cout << std::endl;
367 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
368 void amg_coarse_rs0(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
370 unsigned int total_points;
373 Slicing.slice(level, A, Pointvector);
377 #ifdef VIENNACL_WITH_OPENMP
378 #pragma omp parallel for
380 for (
long i=0; i<static_cast<long>(Slicing.threads_); ++i)
386 Slicing.Offset[i+1][level+1] = Slicing.Pointvector_slice[i][level].get_cpoints();
387 #ifdef VIENNACL_WITH_OPENMP
390 total_points += Slicing.Pointvector_slice[i][level].get_cpoints();
394 if (total_points != 0)
396 #ifdef VIENNACL_WITH_OPENMP
397 #pragma omp parallel for
399 for (
long i=0; i<static_cast<long>(Slicing.threads_); ++i)
402 if (Slicing.Offset[i+1][level+1] == 0)
405 for (
unsigned int j=0; j<Slicing.A_slice[i][level].size1(); ++j)
406 Slicing.Pointvector_slice[i][level].make_cpoint(Slicing.Pointvector_slice[i][level][j]);
407 Slicing.Offset[i+1][level+1] =
static_cast<unsigned int>(Slicing.A_slice[i][level].size1());
412 for (
unsigned int i=2; i<=Slicing.threads_; ++i)
413 Slicing.Offset[i][level+1] += Slicing.Offset[i-1][level+1];
416 Slicing.join(level, Pointvector);
422 #if defined(VIENNACL_AMG_DEBUG)// or defined (VIENNACL_AMG_DEBUGBENCH)
423 for (
unsigned int i=0; i<Slicing._threads; ++i)
425 unsigned int c_points = Slicing.Pointvector_slice[i][level].get_cpoints();
426 unsigned int f_points = Slicing.Pointvector_slice[i][level].get_fpoints();
427 std::cout <<
"Thread " << i <<
": ";
428 std::cout <<
"No of C points = " << c_points <<
", ";
429 std::cout <<
"No of F points = " << f_points << std::endl;
441 template <
typename InternalType1,
typename InternalType2,
typename InternalType3>
442 void amg_coarse_rs3(
unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing,
amg_tag & tag)
452 boost::numeric::ublas::vector<unsigned int> Offset = boost::numeric::ublas::vector<unsigned int> (Slicing.Offset.size());
453 for (i=0; i<Slicing.Offset.size(); ++i)
454 Offset[i] = Slicing.Offset[i][level];
457 for (i=0; i<Slicing.threads_; ++i)
460 for (j=Offset[i]; j<Offset[i+1]; ++j)
462 point1 = Pointvector[level][j];
486 if ((*iter2)->get_index() == (*iter3)->get_index())
492 else if ((*iter2)->get_index() < (*iter3)->get_index())
539 Pointvector[level].switch_ftoc(point2);
541 for (
unsigned int j=i+1; j<=Slicing.threads_; ++j)
542 Slicing.Offset[j][level+1]++;
550 #ifdef VIENNACL_AMG_DEBUG
551 std::cout <<
"After 3rd pass:" << std::endl;
552 std::cout <<
"Coarse Points:" << std::endl;
553 boost::numeric::ublas::vector<bool> C;
554 Pointvector[level].get_C(C);
556 std::cout <<
"Fine Points:" << std::endl;
557 boost::numeric::ublas::vector<bool> F;
558 Pointvector[level].get_F(F);
570 template <
typename InternalType1,
typename InternalType2>
573 typedef typename InternalType1::value_type SparseMatrixType;
574 typedef typename InternalType2::value_type PointVectorType;
575 typedef typename SparseMatrixType::value_type ScalarType;
577 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
578 typedef typename SparseMatrixType::iterator2 InternalColIterator;
585 if (A[level].
size1() == 1)
return;
589 #ifdef VIENNACL_WITH_OPENMP
590 #pragma omp parallel for private (x,y,diag)
592 for (x=0; x<static_cast<long>(A[level].size1()); ++x)
594 InternalRowIterator row_iter = A[level].begin1();
596 diag = A[level](x,x);
597 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
599 y =
static_cast<long>(col_iter.index2());
600 if (y == x || (std::fabs(*col_iter) >= tag.
get_threshold()*pow(0.5, static_cast<double>(level-1)) * std::sqrt(std::fabs(diag*A[level](y,y)))))
608 #ifdef VIENNACL_AMG_DEBUG
609 std::cout <<
"Neighborhoods:" << std::endl;
610 boost::numeric::ublas::matrix<bool> mat;
611 Pointvector[level].get_influence_matrix(mat);
616 for (
typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
639 #ifdef VIENNACL_AMG_DEBUG
640 std::cout <<
"After aggregation:" << std::endl;
641 std::cout <<
"Coarse Points:" << std::endl;
642 boost::numeric::ublas::vector<bool> C;
643 Pointvector[level].get_C(C);
645 std::cout <<
"Fine Points:" << std::endl;
646 boost::numeric::ublas::vector<bool> F;
647 Pointvector[level].get_F(F);
649 std::cout <<
"Aggregates:" << std::endl;
Debug functionality for AMG. To be removed.
void add_influencing_point(amg_point *point)
Definition: amg_base.hpp:867
A class for the AMG points. Saves point index and influence measure Holds information whether point i...
Definition: amg_base.hpp:816
#define VIENNACL_AMG_COARSE_RS
Definition: amg_base.hpp:41
void amg_coarse_classic(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, amg_tag &tag)
Classical (RS) two-pass coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_CLASSIC) ...
Definition: amg_coarse.hpp:245
void amg_coarse_ag(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, amg_tag &tag)
AG (aggregation based) coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_SA)
Definition: amg_coarse.hpp:571
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
#define VIENNACL_AMG_COARSE_RS3
Definition: amg_base.hpp:44
unsigned int get_coarse() const
Definition: amg_base.hpp:90
bool is_undecided() const
Definition: amg_base.hpp:860
void printmatrix(MatrixType &, int)
Definition: amg_debug.hpp:77
void set_aggregate(unsigned int aggregate)
Definition: amg_base.hpp:855
#define VIENNACL_AMG_COARSE_RS0
Definition: amg_base.hpp:43
#define VIENNACL_AMG_COARSE_AG
Definition: amg_base.hpp:45
void make_cpoint()
Definition: amg_base.hpp:890
void amg_coarse_classic_onepass(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, amg_tag &tag)
Classical (RS) one-pass coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_CLASSIC_ONEPASS) ...
Definition: amg_coarse.hpp:157
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:61
iterator begin_influencing()
Definition: amg_base.hpp:907
bool is_cpoint() const
Definition: amg_base.hpp:858
iterator begin_influenced()
Definition: amg_base.hpp:911
void amg_influence(unsigned int level, InternalType1 const &A, InternalType2 &Pointvector, amg_tag &tag)
Determines strong influences in system matrix, classical approach (RS). Multithreaded! ...
Definition: amg_coarse.hpp:71
void make_fpoint()
Definition: amg_base.hpp:897
double get_threshold() const
Definition: amg_base.hpp:96
void printvector(VectorType const &)
Definition: amg_debug.hpp:80
vector_expression< const matrix_base< NumericT, F >, const int, op_matrix_diag > diag(const matrix_base< NumericT, F > &A, int k=0)
Definition: matrix.hpp:895
unsigned int get_index() const
Definition: amg_base.hpp:853
#define VIENNACL_AMG_COARSE_ONEPASS
Definition: amg_base.hpp:42
bool is_influencing(amg_point *point) const
Definition: amg_base.hpp:865
void amg_coarse_rs3(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, InternalType3 &Slicing, amg_tag &tag)
RS3 coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_RS3)
Definition: amg_coarse.hpp:442
Defines an iterator for the sparse vector type.
Definition: amg_base.hpp:203
iterator end_influenced()
Definition: amg_base.hpp:912
Helper classes and functions for the AMG preconditioner. Experimental.
void amg_coarse(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, InternalType3 &Slicing, amg_tag &tag)
Calls the right coarsening procedure.
Definition: amg_coarse.hpp:52
bool is_fpoint() const
Definition: amg_base.hpp:859
iterator end_influencing()
Definition: amg_base.hpp:908
void amg_coarse_rs0(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, InternalType3 &Slicing, amg_tag &tag)
Parallel classical RS0 coarsening. Multi-Threaded! (VIENNACL_AMG_COARSE_RS0 || VIENNACL_AMG_COARSE_RS...
Definition: amg_coarse.hpp:368