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 InternalT1,
typename InternalT2,
typename InternalT3>
52 void amg_coarse(
unsigned int level, InternalT1 & A, InternalT2 & pointvector, InternalT3 & slicing,
amg_tag & tag)
70 template<
typename InternalT1,
typename InternalT2>
73 typedef typename InternalT1::value_type SparseMatrixType;
74 typedef typename InternalT2::value_type PointVectorType;
75 typedef typename SparseMatrixType::value_type
ScalarType;
76 typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
77 typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
79 #ifdef VIENNACL_WITH_OPENMP
80 #pragma omp parallel for
82 for (
long i=0; i<static_cast<long>(A[level].size1()); ++i)
85 if (A[level](static_cast<unsigned int>(i),
static_cast<unsigned int>(i)) < 0)
88 ConstRowIterator row_iter = A[level].begin1();
89 row_iter +=
static_cast<unsigned int>(i);
92 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
94 if (i == (
unsigned int) col_iter.index2())
continue;
96 if (max > *col_iter) max = *col_iter;
98 if (max < *col_iter) max = *col_iter;
102 if (std::fabs(max) <= 0)
106 for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
108 unsigned int j =
static_cast<unsigned int>(col_iter.index2());
116 pointvector[level][
static_cast<unsigned int>(i)]->add_influencing_point(pointvector[level][static_cast<unsigned int>(j)]);
121 #ifdef VIENNACL_AMG_DEBUG
122 std::cout <<
"Influence Matrix: " << std::endl;
123 boost::numeric::ublas::matrix<bool> mat;
124 Pointvector[level].get_influence_matrix(mat);
129 for (
typename PointVectorType::iterator iter = pointvector[level].begin(); iter != pointvector[level].end(); ++iter)
130 for (
typename amg_point::iterator iter2 = (*iter)->begin_influencing(); iter2 != (*iter)->end_influencing(); ++iter2)
131 (*iter2)->add_influenced_point(*iter);
133 #ifdef VIENNACL_AMG_DEBUG
134 std::cout <<
"Influence Measures: " << std::endl;
135 boost::numeric::ublas::vector<unsigned int> temp;
136 pointvector[level].get_influence(temp);
138 std::cout <<
"Point Sorting: " << std::endl;
139 Pointvector[level].get_sorting(temp);
152 template<
typename InternalT1,
typename InternalT2>
161 #ifdef VIENNACL_WITH_OPENMP
162 #pragma omp parallel for
164 for (
long i=0; i<static_cast<long>(pointvector[level].size()); ++i)
165 pointvector[level][static_cast<unsigned int>(i)]->calc_influence();
168 pointvector[level].sort();
171 while ((c_point = pointvector[level].get_nextpoint()) != NULL)
174 pointvector[level].make_cpoint(c_point);
191 pointvector[level].add_influence(point2,1);
213 #if defined (VIENNACL_AMG_DEBUG)// or defined (VIENNACL_AMG_DEBUGBENCH)
214 unsigned int c_points = pointvector[level].get_cpoints();
215 unsigned int f_points = pointvector[level].get_fpoints();
216 std::cout <<
"1st pass: Level " << level <<
": ";
217 std::cout <<
"No of C points = " << c_points <<
", ";
218 std::cout <<
"No of F points = " << f_points << std::endl;
221 #ifdef VIENNACL_AMG_DEBUG
222 std::cout <<
"Coarse Points:" << std::endl;
223 boost::numeric::ublas::vector<bool> C;
224 pointvector[level].get_C(C);
226 std::cout <<
"Fine Points:" << std::endl;
227 boost::numeric::ublas::vector<bool> F;
228 pointvector[level].get_F(F);
239 template<
typename InternalT1,
typename InternalT2>
242 typedef typename InternalT2::value_type PointVectorType;
251 for (
typename PointVectorType::iterator iter = pointvector[level].begin(); iter != pointvector[level].end(); ++iter)
277 if ((*iter2)->get_index() == (*iter3)->get_index())
283 else if ((*iter2)->get_index() < (*iter3)->get_index())
329 #ifdef VIENNACL_AMG_DEBUG
330 std::cout <<
"After 2nd pass:" << std::endl;
331 std::cout <<
"Coarse Points:" << std::endl;
332 boost::numeric::ublas::vector<bool> C;
333 pointvector[level].get_C(C);
335 std::cout <<
"Fine Points:" << std::endl;
336 boost::numeric::ublas::vector<bool> F;
337 pointvector[level].get_F(F);
341 #ifdef VIENNACL_AMG_DEBUG
342 #ifdef VIENNACL_WITH_OPENMP
346 std::cout <<
"No C and no F point: ";
347 for (
typename PointVectorType::iterator iter = pointvector[level].begin(); iter != pointvector[level].end(); ++iter)
348 if ((*iter)->is_undecided())
349 std::cout << (*iter)->get_index() <<
" ";
350 std::cout << std::endl;
362 template<
typename InternalT1,
typename InternalT2,
typename InternalT3>
363 void amg_coarse_rs0(
unsigned int level, InternalT1 & A, InternalT2 & pointvector, InternalT3 & slicing,
amg_tag & tag)
365 unsigned int total_points;
368 slicing.slice(level, A, pointvector);
372 #ifdef VIENNACL_WITH_OPENMP
373 #pragma omp parallel for
375 for (
long i2=0; i2<static_cast<long>(slicing.threads_); ++i2)
377 std::size_t i =
static_cast<std::size_t
>(i2);
382 slicing.offset_[i+1][level+1] = slicing.pointvector_slice_[i][level].get_cpoints();
383 #ifdef VIENNACL_WITH_OPENMP
386 total_points += slicing.pointvector_slice_[i][level].get_cpoints();
390 if (total_points != 0)
392 #ifdef VIENNACL_WITH_OPENMP
393 #pragma omp parallel for
395 for (
long i2=0; i2<static_cast<long>(slicing.threads_); ++i2)
397 std::size_t i =
static_cast<std::size_t
>(i2);
400 if (slicing.offset_[i+1][level+1] == 0)
403 for (
unsigned int j=0; j<slicing.A_slice_[i][level].size1(); ++j)
404 slicing.pointvector_slice_[i][level].make_cpoint(slicing.pointvector_slice_[i][level][j]);
405 slicing.offset_[i+1][level+1] =
static_cast<unsigned int>(slicing.A_slice_[i][level].size1());
410 for (
unsigned int i2=2; i2<=slicing.threads_; ++i2)
412 std::size_t i =
static_cast<std::size_t
>(i2);
413 slicing.offset_[i][level+1] += slicing.offset_[i-1][level+1];
417 slicing.join(level, pointvector);
423 #if defined(VIENNACL_AMG_DEBUG)// or defined (VIENNACL_AMG_DEBUGBENCH)
424 for (
unsigned int i=0; i<slicing._threads; ++i)
426 unsigned int c_points = slicing.pointvector_slice_[i][level].get_cpoints();
427 unsigned int f_points = slicing.pointvector_slice_[i][level].get_fpoints();
428 std::cout <<
"Thread " << i <<
": ";
429 std::cout <<
"No of C points = " << c_points <<
", ";
430 std::cout <<
"No of F points = " << f_points << std::endl;
442 template<
typename InternalT1,
typename InternalT2,
typename InternalT3>
443 void amg_coarse_rs3(
unsigned int level, InternalT1 & A, InternalT2 & pointvector, InternalT3 & slicing,
amg_tag & tag)
453 boost::numeric::ublas::vector<unsigned int> offset = boost::numeric::ublas::vector<unsigned int> (slicing.offset_.size());
454 for (i=0; i<slicing.offset_.size(); ++i)
455 offset[i] = slicing.offset_[i][level];
458 for (i=0; i<slicing.threads_; ++i)
461 for (j=offset[i]; j<offset[i+1]; ++j)
463 point1 = pointvector[level][j];
487 if ((*iter2)->get_index() == (*iter3)->get_index())
493 else if ((*iter2)->get_index() < (*iter3)->get_index())
542 for (
unsigned int k=i+1; k<=slicing.threads_; ++k)
543 slicing.offset_[k][level+1]++;
551 #ifdef VIENNACL_AMG_DEBUG
552 std::cout <<
"After 3rd pass:" << std::endl;
553 std::cout <<
"Coarse Points:" << std::endl;
554 boost::numeric::ublas::vector<bool> C;
555 pointvector[level].get_C(C);
557 std::cout <<
"Fine Points:" << std::endl;
558 boost::numeric::ublas::vector<bool> F;
559 pointvector[level].get_F(F);
571 template<
typename InternalT1,
typename InternalT2>
574 typedef typename InternalT1::value_type SparseMatrixType;
575 typedef typename InternalT2::value_type PointVectorType;
576 typedef typename SparseMatrixType::value_type
ScalarType;
578 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
579 typedef typename SparseMatrixType::iterator2 InternalColIterator;
584 if (A[level].
size1() == 1)
589 #ifdef VIENNACL_WITH_OPENMP
590 #pragma omp parallel for
592 for (
long x=0; x<static_cast<long>(A[level].size1()); ++x)
594 InternalRowIterator row_iter = A[level].begin1();
595 row_iter += std::size_t(x);
596 ScalarType
diag = A[level](
static_cast<unsigned int>(x),static_cast<unsigned int>(x));
597 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
599 long 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](static_cast<unsigned int>(y),static_cast<unsigned int>(y))))))
603 pointvector[level][
static_cast<unsigned int>(x)]->add_influencing_point(pointvector[level][static_cast<unsigned int>(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;
void amg_coarse_rs3(unsigned int level, InternalT1 &A, InternalT2 &pointvector, InternalT3 &slicing, amg_tag &tag)
RS3 coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_RS3)
Debug functionality for AMG. To be removed.
void set_aggregate(unsigned int aggregate)
double get_threshold() const
void amg_coarse_classic(unsigned int level, InternalT1 &A, InternalT2 &pointvector, amg_tag &tag)
Classical (RS) two-pass coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_CLASSIC) ...
A class for the AMG points. Saves point index and influence measure Holds information whether point i...
#define VIENNACL_AMG_COARSE_RS
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
#define VIENNACL_AMG_COARSE_RS3
void amg_coarse_classic_onepass(unsigned int level, InternalT1 &A, InternalT2 &pointvector, amg_tag &tag)
Classical (RS) one-pass coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_CLASSIC_ONEPASS) ...
T max(const T &lhs, const T &rhs)
Maximum.
#define VIENNACL_AMG_COARSE_RS0
#define VIENNACL_AMG_COARSE_AG
iterator begin_influenced()
bool is_influencing(amg_point *point) const
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
iterator begin_influencing()
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
void printvector(VectorT const &)
bool is_undecided() const
void amg_coarse(unsigned int level, InternalT1 &A, InternalT2 &pointvector, InternalT3 &slicing, amg_tag &tag)
Calls the right coarsening procedure.
unsigned int get_coarse() const
void amg_coarse_ag(unsigned int level, InternalT1 &A, InternalT2 &pointvector, amg_tag &tag)
AG (aggregation based) coarsening. Single-Threaded! (VIENNACL_AMG_COARSE_SA)
#define VIENNACL_AMG_COARSE_ONEPASS
void amg_coarse_rs0(unsigned int level, InternalT1 &A, InternalT2 &pointvector, InternalT3 &slicing, amg_tag &tag)
Parallel classical RS0 coarsening. Multi-Threaded! (VIENNACL_AMG_COARSE_RS0 || VIENNACL_AMG_COARSE_RS...
Defines an iterator for the sparse vector type.
iterator end_influenced()
unsigned int get_index() const
Helper classes and functions for the AMG preconditioner. Experimental.
iterator end_influencing()
void printmatrix(MatrixT &, int)
void amg_influence(unsigned int level, InternalT1 const &A, InternalT2 &pointvector, amg_tag &tag)
Determines strong influences in system matrix, classical approach (RS). Multithreaded! ...