1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
25 #include <boost/numeric/ublas/vector.hpp>
30 #ifdef VIENNACL_WITH_OPENMP
52 template<
typename InternalT1,
typename InternalT2>
53 void amg_interpol(
unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector,
amg_tag & tag)
71 template<
typename InternalT1,
typename InternalT2>
74 typedef typename InternalT1::value_type SparseMatrixType;
76 typedef typename SparseMatrixType::value_type
ScalarType;
77 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
78 typedef typename SparseMatrixType::iterator2 InternalColIterator;
80 unsigned int c_points = pointvector[level].get_cpoints();
83 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
87 pointvector[level].build_index();
90 #ifdef VIENNACL_WITH_OPENMP
91 #pragma omp parallel for
93 for (
long x=0; x < static_cast<long>(pointvector[level].size()); ++x)
95 amg_point *pointx = pointvector[level][
static_cast<unsigned int>(x)];
109 InternalRowIterator row_iter = A[level].begin1();
113 ScalarType row_sum = 0;
114 ScalarType c_sum = 0;
116 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
118 long y =
static_cast<long>(col_iter.index2());
126 row_sum += *col_iter;
128 amg_point *pointy = pointvector[level][
static_cast<unsigned int>(y)];
134 ScalarType temp_res = -row_sum/(c_sum*
diag);
143 if (temp_res > 0 || temp_res < 0)
144 P[level](
static_cast<unsigned int>(x), pointy->
get_coarse_index()) = temp_res * A[level](static_cast<unsigned int>(x),pointy->
get_index());
157 #ifdef VIENNACL_AMG_DEBUG
158 std::cout <<
"Prolongation Matrix:" << std::endl;
170 template<
typename InternalT1,
typename InternalT2>
173 typedef typename InternalT1::value_type SparseMatrixType;
175 typedef typename SparseMatrixType::value_type
ScalarType;
176 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
177 typedef typename SparseMatrixType::iterator2 InternalColIterator;
179 unsigned int c_points = pointvector[level].get_cpoints();
182 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
186 pointvector[level].build_index();
189 #ifdef VIENNACL_WITH_OPENMP
190 #pragma omp parallel for
192 for (
long x=0; x < static_cast<long>(pointvector[level].size()); ++x)
194 amg_point *pointx = pointvector[level][
static_cast<unsigned int>(x)];
195 int diag_sign = (A[level](
static_cast<unsigned int>(x),static_cast<unsigned int>(x)) > 0) ? 1 : -1;
205 InternalRowIterator row_iter = A[level].begin1();
208 ScalarType weak_sum = 0;
211 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
213 long k =
static_cast<unsigned int>(col_iter.index2());
214 amg_point *pointk = pointvector[level][
static_cast<unsigned int>(k)];
219 weak_sum += *col_iter;
233 if (A[level](static_cast<unsigned int>(k),static_cast<unsigned int>(m)) *
ScalarType(diag_sign) < 0)
234 c_sum_row[static_cast<unsigned int>(k)] += A[level](
static_cast<unsigned int>(k), static_cast<unsigned int>(m));
249 ScalarType strong_sum = 0;
253 long k = iter2.index();
255 if (A[level](static_cast<unsigned int>(k),static_cast<unsigned int>(y)) *
ScalarType(diag_sign) < 0)
256 strong_sum += (A[level](
static_cast<unsigned int>(x),static_cast<unsigned int>(k)) * A[level](static_cast<unsigned int>(k),
static_cast<unsigned int>(y))) / (*iter2);
260 ScalarType temp_res = - (A[level](
static_cast<unsigned int>(x),static_cast<unsigned int>(y)) + strong_sum) / (weak_sum);
261 if (temp_res < 0 || temp_res > 0)
262 P[level](
static_cast<unsigned int>(x),pointy->
get_coarse_index()) = temp_res;
272 #ifdef VIENNACL_AMG_DEBUG
273 std::cout <<
"Prolongation Matrix:" << std::endl;
284 template<
typename SparseMatrixT>
287 typedef typename SparseMatrixT::value_type
ScalarType;
288 typedef typename SparseMatrixT::iterator1 InternalRowIterator;
289 typedef typename SparseMatrixT::iterator2 InternalColIterator;
291 ScalarType row_max, row_min, row_sum_pos, row_sum_neg, row_sum_pos_scale, row_sum_neg_scale;
293 InternalRowIterator row_iter = P.begin1();
303 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
305 if (*col_iter > row_max)
307 if (*col_iter < row_min)
310 row_sum_pos += *col_iter;
312 row_sum_neg += *col_iter;
315 row_sum_pos_scale = row_sum_pos;
316 row_sum_neg_scale = row_sum_neg;
319 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
323 row_sum_pos_scale -= *col_iter;
328 row_sum_pos_scale -= *col_iter;
334 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
337 *col_iter = *col_iter *(row_sum_pos/row_sum_pos_scale);
339 *col_iter = *col_iter *(row_sum_neg/row_sum_neg_scale);
349 template<
typename InternalT1,
typename InternalT2>
352 typedef typename InternalT1::value_type SparseMatrixType;
358 unsigned int c_points = pointvector[level].get_cpoints();
360 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
364 pointvector[level].build_index();
367 #ifdef VIENNACL_WITH_OPENMP
368 #pragma omp parallel for
370 for (
long x=0; x<static_cast<long>(pointvector[level].size()); ++x)
372 amg_point *pointx = pointvector[level][
static_cast<unsigned int>(x)];
375 P[level](
static_cast<unsigned int>(x), pointy->get_coarse_index()) = 1;
378 #ifdef VIENNACL_AMG_DEBUG
379 std::cout <<
"Aggregation based Prolongation:" << std::endl;
391 template<
typename InternalT1,
typename InternalT2>
394 typedef typename InternalT1::value_type SparseMatrixType;
396 typedef typename SparseMatrixType::value_type
ScalarType;
397 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
398 typedef typename SparseMatrixType::iterator2 InternalColIterator;
400 unsigned int c_points = pointvector[level].get_cpoints();
402 InternalT1 P_tentative = InternalT1(P.size());
403 SparseMatrixType Jacobi = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), static_cast<unsigned int>(A[level].
size2()));
405 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
409 #ifdef VIENNACL_WITH_OPENMP
410 #pragma omp parallel for
412 for (
long x=0; x<static_cast<long>(A[level].size1()); ++x)
415 InternalRowIterator row_iter = A[level].begin1();
417 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
419 long y =
static_cast<long>(col_iter.index2());
426 else if (!pointvector[level][static_cast<unsigned int>(x)]->is_influencing(pointvector[level][static_cast<unsigned int>(y)]))
429 Jacobi (static_cast<unsigned int>(x), static_cast<unsigned int>(y)) = *col_iter;
431 InternalRowIterator row_iter2 = Jacobi.begin1();
434 for (InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
439 Jacobi (static_cast<unsigned int>(x), static_cast<unsigned int>(x)) = 1 -
static_cast<ScalarType
>(tag.
get_interpolweight());
442 #ifdef VIENNACL_AMG_DEBUG
443 std::cout <<
"Jacobi Matrix:" << std::endl;
450 #ifdef VIENNACL_AMG_DEBUG
451 std::cout <<
"Tentative Prolongation:" << std::endl;
458 #ifdef VIENNACL_AMG_DEBUG
459 std::cout <<
"Prolongation Matrix:" << std::endl;
#define VIENNACL_AMG_INTERPOL_SA
Debug functionality for AMG. To be removed.
#define VIENNACL_AMG_INTERPOL_AG
#define VIENNACL_AMG_INTERPOL_DIRECT
A class for the AMG points. Saves point index and influence measure Holds information whether point i...
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
unsigned int get_aggregate()
A class for the sparse vector type.
void amg_interpol_ag(unsigned int level, InternalT1 &A, InternalT1 &P, InternalT2 &pointvector, amg_tag)
AG (aggregation based) interpolation. Multi-Threaded! (VIENNACL_INTERPOL_SA)
unsigned int get_coarse_index() const
void amg_interpol_sa(unsigned int level, InternalT1 &A, InternalT1 &P, InternalT2 &pointvector, amg_tag &tag)
SA (smoothed aggregate) interpolation. Multi-Threaded! (VIENNACL_INTERPOL_SA)
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
bool is_influencing(amg_point *point) const
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
void amg_interpol(unsigned int level, InternalT1 &A, InternalT1 &P, InternalT2 &pointvector, amg_tag &tag)
Calls the right function to build interpolation matrix.
unsigned int get_interpol() const
iterator begin_influencing()
void amg_truncate_row(SparseMatrixT &P, unsigned int row, amg_tag &tag)
Interpolation truncation (for VIENNACL_AMG_INTERPOL_DIRECT and VIENNACL_AMG_INTERPOL_CLASSIC) ...
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
double get_interpolweight() const
#define VIENNACL_AMG_INTERPOL_CLASSIC
Defines an iterator for the sparse vector type.
unsigned int get_index() const
void amg_mat_prod(SparseMatrixT &A, SparseMatrixT &B, SparseMatrixT &RES)
Sparse matrix product. Calculates RES = A*B.
Helper classes and functions for the AMG preconditioner. Experimental.
iterator end_influencing()
void printmatrix(MatrixT &, int)
void amg_interpol_direct(unsigned int level, InternalT1 &A, InternalT1 &P, InternalT2 &pointvector, amg_tag &tag)
Direct interpolation. Multi-threaded! (VIENNACL_AMG_INTERPOL_DIRECT)
void amg_interpol_classic(unsigned int level, InternalT1 &A, InternalT1 &P, InternalT2 &pointvector, amg_tag &tag)
Classical interpolation. Don't use with onepass classical coarsening or RS0 (Yang, p.14)! Multi-threaded! (VIENNACL_AMG_INTERPOL_CLASSIC)