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 InternalType1,
typename InternalType2>
53 void amg_interpol(
unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
70 template <
typename InternalType1,
typename InternalType2>
73 typedef typename InternalType1::value_type SparseMatrixType;
75 typedef typename SparseMatrixType::value_type ScalarType;
76 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
77 typedef typename SparseMatrixType::iterator2 InternalColIterator;
80 ScalarType row_sum, c_sum,
diag;
84 unsigned int c_points = Pointvector[level].get_cpoints();
87 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()),c_points);
91 Pointvector[level].build_index();
94 #ifdef VIENNACL_WITH_OPENMP
95 #pragma omp parallel for private (pointx,pointy,row_sum,c_sum,temp_res,y,x,diag)
97 for (x=0; x < static_cast<long>(Pointvector[level].size()); ++x)
99 pointx = Pointvector[level][x];
113 InternalRowIterator row_iter = A[level].begin1();
117 row_sum = c_sum = diag = 0;
118 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
120 y =
static_cast<long>(col_iter.index2());
128 row_sum += *col_iter;
130 pointy = Pointvector[level][y];
136 temp_res = -row_sum/(c_sum*
diag);
159 #ifdef VIENNACL_AMG_DEBUG
160 std::cout <<
"Prolongation Matrix:" << std::endl;
172 template <
typename InternalType1,
typename InternalType2>
175 typedef typename InternalType1::value_type SparseMatrixType;
177 typedef typename SparseMatrixType::value_type ScalarType;
178 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
179 typedef typename SparseMatrixType::iterator2 InternalColIterator;
182 ScalarType weak_sum, strong_sum;
185 amg_point *pointx, *pointy, *pointk, *pointm;
188 unsigned int c_points = Pointvector[level].get_cpoints();
191 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
195 Pointvector[level].build_index();
198 #ifdef VIENNACL_WITH_OPENMP
199 #pragma omp parallel for private (pointx,pointy,pointk,pointm,weak_sum,strong_sum,c_sum_row,temp_res,x,y,k,m,diag_sign)
201 for (x=0; x < static_cast<long>(Pointvector[level].size()); ++x)
203 pointx = Pointvector[level][x];
204 if (A[level](x,x) > 0)
217 InternalRowIterator row_iter = A[level].begin1();
223 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
225 k =
static_cast<unsigned int>(col_iter.index2());
226 pointk = Pointvector[level][k];
231 weak_sum += *col_iter;
245 if (A[level](k,m) * diag_sign < 0)
246 c_sum_row[k] += A[level](k,m);
267 if (A[level](k,y) * diag_sign < 0)
268 strong_sum += (A[level](x,k) * A[level](k,y)) / (*iter2);
272 temp_res = - (A[level](x,y) + strong_sum) / (weak_sum);
284 #ifdef VIENNACL_AMG_DEBUG
285 std::cout <<
"Prolongation Matrix:" << std::endl;
296 template <
typename SparseMatrixType>
299 typedef typename SparseMatrixType::value_type ScalarType;
300 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
301 typedef typename SparseMatrixType::iterator2 InternalColIterator;
303 ScalarType row_max, row_min, row_sum_pos, row_sum_neg, row_sum_pos_scale, row_sum_neg_scale;
305 InternalRowIterator row_iter = P.begin1();
315 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
317 if (*col_iter > row_max)
319 if (*col_iter < row_min)
322 row_sum_pos += *col_iter;
324 row_sum_neg += *col_iter;
327 row_sum_pos_scale = row_sum_pos;
328 row_sum_neg_scale = row_sum_neg;
331 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
335 row_sum_pos_scale -= *col_iter;
340 row_sum_pos_scale -= *col_iter;
346 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
349 *col_iter = *col_iter *(row_sum_pos/row_sum_pos_scale);
351 *col_iter = *col_iter *(row_sum_neg/row_sum_neg_scale);
361 template <
typename InternalType1,
typename InternalType2>
364 typedef typename InternalType1::value_type SparseMatrixType;
372 unsigned int c_points = Pointvector[level].get_cpoints();
374 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
378 Pointvector[level].build_index();
381 #ifdef VIENNACL_WITH_OPENMP
382 #pragma omp parallel for private (x,pointx)
384 for (x=0; x<static_cast<long>(Pointvector[level].size()); ++x)
386 pointx = Pointvector[level][x];
392 #ifdef VIENNACL_AMG_DEBUG
393 std::cout <<
"Aggregation based Prolongation:" << std::endl;
405 template <
typename InternalType1,
typename InternalType2>
406 void amg_interpol_sa(
unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector,
amg_tag & tag)
408 typedef typename InternalType1::value_type SparseMatrixType;
410 typedef typename SparseMatrixType::value_type ScalarType;
411 typedef typename SparseMatrixType::iterator1 InternalRowIterator;
412 typedef typename SparseMatrixType::iterator2 InternalColIterator;
416 unsigned int c_points = Pointvector[level].get_cpoints();
418 InternalType1 P_tentative = InternalType1(P.size());
419 SparseMatrixType Jacobi = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), static_cast<unsigned int>(A[level].
size2()));
421 P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].
size1()), c_points);
425 #ifdef VIENNACL_WITH_OPENMP
426 #pragma omp parallel for private (x,y,diag)
428 for (x=0; x<static_cast<long>(A[level].size1()); ++x)
431 InternalRowIterator row_iter = A[level].begin1();
433 for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
435 y =
static_cast<long>(col_iter.index2());
442 else if (!Pointvector[level][x]->is_influencing(Pointvector[level][y]))
445 Jacobi (x,y) = *col_iter;
447 InternalRowIterator row_iter2 = Jacobi.begin1();
450 for (InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
458 #ifdef VIENNACL_AMG_DEBUG
459 std::cout <<
"Jacobi Matrix:" << std::endl;
466 #ifdef VIENNACL_AMG_DEBUG
467 std::cout <<
"Tentative Prolongation:" << std::endl;
474 #ifdef VIENNACL_AMG_DEBUG
475 std::cout <<
"Prolongation Matrix:" << std::endl;
void amg_interpol(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Calls the right function to build interpolation matrix.
Definition: amg_interpol.hpp:53
#define VIENNACL_AMG_INTERPOL_SA
Definition: amg_base.hpp:49
iterator begin()
Definition: amg_base.hpp:352
Debug functionality for AMG. To be removed.
void amg_interpol_direct(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Direct interpolation. Multi-threaded! (VIENNACL_AMG_INTERPOL_DIRECT)
Definition: amg_interpol.hpp:71
iterator end()
Definition: amg_base.hpp:354
#define VIENNACL_AMG_INTERPOL_AG
Definition: amg_base.hpp:48
#define VIENNACL_AMG_INTERPOL_DIRECT
Definition: amg_base.hpp:46
A class for the AMG points. Saves point index and influence measure Holds information whether point i...
Definition: amg_base.hpp:816
void amg_interpol_classic(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Classical interpolation. Don't use with onepass classical coarsening or RS0 (Yang, p.14)! Multi-threaded! (VIENNACL_AMG_INTERPOL_CLASSIC)
Definition: amg_interpol.hpp:173
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
A class for the sparse vector type.
Definition: amg_base.hpp:254
unsigned int get_coarse_index() const
Definition: amg_base.hpp:877
void printmatrix(MatrixType &, int)
Definition: amg_debug.hpp:77
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:245
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
unsigned int get_aggregate()
Definition: amg_base.hpp:856
void amg_interpol_ag(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag)
AG (aggregation based) interpolation. Multi-Threaded! (VIENNACL_INTERPOL_SA)
Definition: amg_interpol.hpp:362
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
void amg_interpol_sa(unsigned int level, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
SA (smoothed aggregate) interpolation. Multi-Threaded! (VIENNACL_INTERPOL_SA)
Definition: amg_interpol.hpp:406
void amg_truncate_row(SparseMatrixType &P, unsigned int row, amg_tag &tag)
Interpolation truncation (for VIENNACL_AMG_INTERPOL_DIRECT and VIENNACL_AMG_INTERPOL_CLASSIC) ...
Definition: amg_interpol.hpp:297
void clear()
Definition: amg_base.hpp:288
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
double get_interpolweight() const
Definition: amg_base.hpp:99
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
unsigned int get_interpol() const
Definition: amg_base.hpp:93
void amg_mat_prod(SparseMatrixType &A, SparseMatrixType &B, SparseMatrixType &RES)
Sparse matrix product. Calculates RES = A*B.
Definition: amg_base.hpp:1319
bool is_influencing(amg_point *point) const
Definition: amg_base.hpp:865
#define VIENNACL_AMG_INTERPOL_CLASSIC
Definition: amg_base.hpp:47
Defines an iterator for the sparse vector type.
Definition: amg_base.hpp:203
Helper classes and functions for the AMG preconditioner. Experimental.
bool is_fpoint() const
Definition: amg_base.hpp:859
iterator end_influencing()
Definition: amg_base.hpp:908