ViennaCL - The Vienna Computing Library  1.6.2
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
amg_interpol.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_INTERPOL_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
25 #include <boost/numeric/ublas/vector.hpp>
26 #include <cmath>
28 
29 #include <map>
30 #ifdef VIENNACL_WITH_OPENMP
31 #include <omp.h>
32 #endif
33 
35 
36 namespace viennacl
37 {
38 namespace linalg
39 {
40 namespace detail
41 {
42 namespace amg
43 {
44 
52 template<typename InternalT1, typename InternalT2>
53 void amg_interpol(unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector, amg_tag & tag)
54 {
55  switch (tag.get_interpol())
56  {
57  case VIENNACL_AMG_INTERPOL_DIRECT: amg_interpol_direct (level, A, P, pointvector, tag); break;
58  case VIENNACL_AMG_INTERPOL_CLASSIC: amg_interpol_classic(level, A, P, pointvector, tag); break;
59  case VIENNACL_AMG_INTERPOL_AG: amg_interpol_ag (level, A, P, pointvector, tag); break;
60  case VIENNACL_AMG_INTERPOL_SA: amg_interpol_sa (level, A, P, pointvector, tag); break;
61  }
62 }
63 
71 template<typename InternalT1, typename InternalT2>
72 void amg_interpol_direct(unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector, amg_tag & tag)
73 {
74  typedef typename InternalT1::value_type SparseMatrixType;
75  //typedef typename InternalType2::value_type PointVectorType;
76  typedef typename SparseMatrixType::value_type ScalarType;
77  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
78  typedef typename SparseMatrixType::iterator2 InternalColIterator;
79 
80  unsigned int c_points = pointvector[level].get_cpoints();
81 
82  // Setup Prolongation/Interpolation matrix
83  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
84  P[level].clear();
85 
86  // Assign indices to C points
87  pointvector[level].build_index();
88 
89  // Direct Interpolation (Yang, p.14)
90 #ifdef VIENNACL_WITH_OPENMP
91  #pragma omp parallel for
92 #endif
93  for (long x=0; x < static_cast<long>(pointvector[level].size()); ++x)
94  {
95  amg_point *pointx = pointvector[level][static_cast<unsigned int>(x)];
96  /*if (A[level](x,x) > 0)
97  diag_sign = 1;
98  else
99  diag_sign = -1;*/
100 
101  // When the current line corresponds to a C point then the diagonal coefficient is 1 and the rest 0
102  if (pointx->is_cpoint())
103  P[level](static_cast<unsigned int>(x),pointx->get_coarse_index()) = 1;
104 
105  // When the current line corresponds to a F point then the diagonal is 0 and the rest has to be computed (Yang, p.14)
106  if (pointx->is_fpoint())
107  {
108  // Jump to row x
109  InternalRowIterator row_iter = A[level].begin1();
110  row_iter += vcl_size_t(x);
111 
112  // Row sum of coefficients (without diagonal) and sum of influencing C point coefficients has to be computed
113  ScalarType row_sum = 0;
114  ScalarType c_sum = 0;
115  ScalarType diag = 0;
116  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
117  {
118  long y = static_cast<long>(col_iter.index2());
119  if (x == y)// || *col_iter * diag_sign > 0)
120  {
121  diag += *col_iter;
122  continue;
123  }
124 
125  // Sum all other coefficients in line x
126  row_sum += *col_iter;
127 
128  amg_point *pointy = pointvector[level][static_cast<unsigned int>(y)];
129  // Sum all coefficients that correspond to a strongly influencing C point
130  if (pointy->is_cpoint())
131  if (pointx->is_influencing(pointy))
132  c_sum += *col_iter;
133  }
134  ScalarType temp_res = -row_sum/(c_sum*diag);
135 
136  // Iterate over all strongly influencing points of point x
137  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
138  {
139  amg_point *pointy = *iter;
140  // The value is only non-zero for columns that correspond to a C point
141  if (pointy->is_cpoint())
142  {
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());
145  }
146  }
147 
148  //Truncate interpolation if chosen
149  if (tag.get_interpolweight() > 0)
150  amg_truncate_row(P[level], static_cast<unsigned int>(x), tag);
151  }
152  }
153 
154  // P test
155  //test_interpolation(A[level], P[level], Pointvector[level]);
156 
157  #ifdef VIENNACL_AMG_DEBUG
158  std::cout << "Prolongation Matrix:" << std::endl;
159  printmatrix (P[level]);
160  #endif
161 }
162 
170 template<typename InternalT1, typename InternalT2>
171 void amg_interpol_classic(unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector, amg_tag & tag)
172 {
173  typedef typename InternalT1::value_type SparseMatrixType;
174  //typedef typename InternalType2::value_type PointVectorType;
175  typedef typename SparseMatrixType::value_type ScalarType;
176  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
177  typedef typename SparseMatrixType::iterator2 InternalColIterator;
178 
179  unsigned int c_points = pointvector[level].get_cpoints();
180 
181  // Setup Prolongation/Interpolation matrix
182  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
183  P[level].clear();
184 
185  // Assign indices to C points
186  pointvector[level].build_index();
187 
188  // Classical Interpolation (Yang, p.13-14)
189 #ifdef VIENNACL_WITH_OPENMP
190  #pragma omp parallel for
191 #endif
192  for (long x=0; x < static_cast<long>(pointvector[level].size()); ++x)
193  {
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;
196 
197  // When the current line corresponds to a C point then the diagonal coefficient is 1 and the rest 0
198  if (pointx->is_cpoint())
199  P[level](static_cast<unsigned int>(x),pointx->get_coarse_index()) = 1;
200 
201  // When the current line corresponds to a F point then the diagonal is 0 and the rest has to be computed (Yang, p.14)
202  if (pointx->is_fpoint())
203  {
204  // Jump to row x
205  InternalRowIterator row_iter = A[level].begin1();
206  row_iter += vcl_size_t(x);
207 
208  ScalarType weak_sum = 0;
209  amg_sparsevector<ScalarType> c_sum_row(static_cast<unsigned int>(A[level].size1()));
210  c_sum_row.clear();
211  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
212  {
213  long k = static_cast<unsigned int>(col_iter.index2());
214  amg_point *pointk = pointvector[level][static_cast<unsigned int>(k)];
215 
216  // Sum of weakly influencing neighbors + diagonal coefficient
217  if (x == k || !pointx->is_influencing(pointk))// || *col_iter * diag_sign > 0)
218  {
219  weak_sum += *col_iter;
220  continue;
221  }
222 
223  // Sums of coefficients in row k (strongly influening F neighbors) of C point neighbors of x are calculated
224  if (pointk->is_fpoint() && pointx->is_influencing(pointk))
225  {
226  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
227  {
228  amg_point *pointm = *iter;
229  long m = pointm->get_index();
230 
231  if (pointm->is_cpoint())
232  // Only use coefficients that have opposite sign of diagonal.
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));
235  }
236  continue;
237  }
238  }
239 
240  // Iterate over all strongly influencing points of point x
241  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
242  {
243  amg_point *pointy = *iter;
244  long y = pointy->get_index();
245 
246  // The value is only non-zero for columns that correspond to a C point
247  if (pointy->is_cpoint())
248  {
249  ScalarType strong_sum = 0;
250  // Calculate term for strongly influencing F neighbors
251  for (typename amg_sparsevector<ScalarType>::iterator iter2 = c_sum_row.begin(); iter2 != c_sum_row.end(); ++iter2)
252  {
253  long k = iter2.index();
254  // Only use coefficients that have opposite sign of diagonal.
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);
257  }
258 
259  // Calculate coefficient
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;
263  }
264  }
265 
266  //Truncate iteration if chosen
267  if (tag.get_interpolweight() > 0)
268  amg_truncate_row(P[level], static_cast<unsigned int>(x), tag);
269  }
270  }
271 
272  #ifdef VIENNACL_AMG_DEBUG
273  std::cout << "Prolongation Matrix:" << std::endl;
274  printmatrix (P[level]);
275  #endif
276 }
277 
284 template<typename SparseMatrixT>
285 void amg_truncate_row(SparseMatrixT & P, unsigned int row, amg_tag & tag)
286 {
287  typedef typename SparseMatrixT::value_type ScalarType;
288  typedef typename SparseMatrixT::iterator1 InternalRowIterator;
289  typedef typename SparseMatrixT::iterator2 InternalColIterator;
290 
291  ScalarType row_max, row_min, row_sum_pos, row_sum_neg, row_sum_pos_scale, row_sum_neg_scale;
292 
293  InternalRowIterator row_iter = P.begin1();
294  row_iter += row;
295 
296  row_max = 0;
297  row_min = 0;
298  row_sum_pos = 0;
299  row_sum_neg = 0;
300 
301  // Truncate interpolation by making values to zero that are a lot smaller than the biggest value in a row
302  // Determine max entry and sum of row (seperately for negative and positive entries)
303  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
304  {
305  if (*col_iter > row_max)
306  row_max = *col_iter;
307  if (*col_iter < row_min)
308  row_min = *col_iter;
309  if (*col_iter > 0)
310  row_sum_pos += *col_iter;
311  if (*col_iter < 0)
312  row_sum_neg += *col_iter;
313  }
314 
315  row_sum_pos_scale = row_sum_pos;
316  row_sum_neg_scale = row_sum_neg;
317 
318  // Make certain values to zero (seperately for negative and positive entries)
319  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
320  {
321  if (*col_iter > 0 && *col_iter < tag.get_interpolweight() * row_max)
322  {
323  row_sum_pos_scale -= *col_iter;
324  *col_iter = 0;
325  }
326  if (*col_iter < 0 && *col_iter > tag.get_interpolweight() * row_min)
327  {
328  row_sum_pos_scale -= *col_iter;
329  *col_iter = 0;
330  }
331  }
332 
333  // Scale remaining values such that row sum is unchanged
334  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
335  {
336  if (*col_iter > 0)
337  *col_iter = *col_iter *(row_sum_pos/row_sum_pos_scale);
338  if (*col_iter < 0)
339  *col_iter = *col_iter *(row_sum_neg/row_sum_neg_scale);
340  }
341 }
342 
349 template<typename InternalT1, typename InternalT2>
350 void amg_interpol_ag(unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector, amg_tag)
351 {
352  typedef typename InternalT1::value_type SparseMatrixType;
353  //typedef typename InternalType2::value_type PointVectorType;
354  //typedef typename SparseMatrixType::value_type ScalarType;
355  //typedef typename SparseMatrixType::iterator1 InternalRowIterator;
356  //typedef typename SparseMatrixType::iterator2 InternalColIterator;
357 
358  unsigned int c_points = pointvector[level].get_cpoints();
359 
360  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
361  P[level].clear();
362 
363  // Assign indices to C points
364  pointvector[level].build_index();
365 
366  // Set prolongation such that F point is interpolated (weight=1) by the aggregate it belongs to (Vanek et al p.6)
367 #ifdef VIENNACL_WITH_OPENMP
368  #pragma omp parallel for
369 #endif
370  for (long x=0; x<static_cast<long>(pointvector[level].size()); ++x)
371  {
372  amg_point *pointx = pointvector[level][static_cast<unsigned int>(x)];
373  amg_point *pointy = pointvector[level][pointx->get_aggregate()];
374  // Point x belongs to aggregate y.
375  P[level](static_cast<unsigned int>(x), pointy->get_coarse_index()) = 1;
376  }
377 
378  #ifdef VIENNACL_AMG_DEBUG
379  std::cout << "Aggregation based Prolongation:" << std::endl;
380  printmatrix(P[level]);
381  #endif
382 }
383 
391 template<typename InternalT1, typename InternalT2>
392 void amg_interpol_sa(unsigned int level, InternalT1 & A, InternalT1 & P, InternalT2 & pointvector, amg_tag & tag)
393 {
394  typedef typename InternalT1::value_type SparseMatrixType;
395  //typedef typename InternalType2::value_type PointVectorType;
396  typedef typename SparseMatrixType::value_type ScalarType;
397  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
398  typedef typename SparseMatrixType::iterator2 InternalColIterator;
399 
400  unsigned int c_points = pointvector[level].get_cpoints();
401 
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()));
404  Jacobi.clear();
405  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
406  P[level].clear();
407 
408  // Build Jacobi Matrix via filtered A matrix (Vanek et al. p.6)
409 #ifdef VIENNACL_WITH_OPENMP
410  #pragma omp parallel for
411 #endif
412  for (long x=0; x<static_cast<long>(A[level].size1()); ++x)
413  {
414  ScalarType diag = 0;
415  InternalRowIterator row_iter = A[level].begin1();
416  row_iter += vcl_size_t(x);
417  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
418  {
419  long y = static_cast<long>(col_iter.index2());
420  // Determine the structure of the Jacobi matrix by using a filtered matrix of A:
421  // The diagonal consists of the diagonal coefficient minus all coefficients of points not in the neighborhood of x.
422  // All other coefficients are the same as in A.
423  // Already use Jacobi matrix to save filtered A matrix to speed up computation.
424  if (x == y)
425  diag += *col_iter;
426  else if (!pointvector[level][static_cast<unsigned int>(x)]->is_influencing(pointvector[level][static_cast<unsigned int>(y)]))
427  diag += -*col_iter;
428  else
429  Jacobi (static_cast<unsigned int>(x), static_cast<unsigned int>(y)) = *col_iter;
430  }
431  InternalRowIterator row_iter2 = Jacobi.begin1();
432  row_iter2 += vcl_size_t(x);
433  // Traverse through filtered A matrix and compute the Jacobi filtering
434  for (InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
435  {
436  *col_iter2 = - static_cast<ScalarType>(tag.get_interpolweight())/diag * *col_iter2;
437  }
438  // Diagonal can be computed seperately.
439  Jacobi (static_cast<unsigned int>(x), static_cast<unsigned int>(x)) = 1 - static_cast<ScalarType>(tag.get_interpolweight());
440  }
441 
442  #ifdef VIENNACL_AMG_DEBUG
443  std::cout << "Jacobi Matrix:" << std::endl;
444  printmatrix(Jacobi);
445  #endif
446 
447  // Use AG interpolation as tentative prolongation
448  amg_interpol_ag(level, A, P_tentative, pointvector, tag);
449 
450  #ifdef VIENNACL_AMG_DEBUG
451  std::cout << "Tentative Prolongation:" << std::endl;
452  printmatrix(P_tentative[level]);
453  #endif
454 
455  // Multiply Jacobi matrix with tentative prolongation to get actual prolongation
456  amg_mat_prod(Jacobi,P_tentative[level],P[level]);
457 
458  #ifdef VIENNACL_AMG_DEBUG
459  std::cout << "Prolongation Matrix:" << std::endl;
460  printmatrix(P[level]);
461  #endif
462 }
463 
464 } //namespace amg
465 } //namespace detail
466 } //namespace linalg
467 } //namespace viennacl
468 
469 #endif
#define VIENNACL_AMG_INTERPOL_SA
Definition: amg_base.hpp:49
Debug functionality for AMG. To be removed.
#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:827
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:260
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
Definition: amg_base.hpp:888
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...
Definition: size.hpp:245
bool is_influencing(amg_point *point) const
Definition: amg_base.hpp:876
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:61
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
Definition: amg_base.hpp:93
void amg_truncate_row(SparseMatrixT &P, unsigned int row, amg_tag &tag)
Interpolation truncation (for VIENNACL_AMG_INTERPOL_DIRECT and VIENNACL_AMG_INTERPOL_CLASSIC) ...
std::size_t vcl_size_t
Definition: forwards.h:74
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
Definition: matrix.hpp:838
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:853
float ScalarType
Definition: fft_1d.cpp:42
#define VIENNACL_AMG_INTERPOL_CLASSIC
Definition: amg_base.hpp:47
Defines an iterator for the sparse vector type.
Definition: amg_base.hpp:209
void amg_mat_prod(SparseMatrixT &A, SparseMatrixT &B, SparseMatrixT &RES)
Sparse matrix product. Calculates RES = A*B.
Definition: amg_base.hpp:1330
Helper classes and functions for the AMG preconditioner. Experimental.
void printmatrix(MatrixT &, int)
Definition: amg_debug.hpp:77
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)