ViennaCL - The Vienna Computing Library  1.5.2
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 InternalType1, typename InternalType2>
53  void amg_interpol(unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & 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  }
70  template <typename InternalType1, typename InternalType2>
71  void amg_interpol_direct(unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
72  {
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::iterator1 InternalRowIterator;
77  typedef typename SparseMatrixType::iterator2 InternalColIterator;
78 
79  ScalarType temp_res;
80  ScalarType row_sum, c_sum, diag;
81  //int diag_sign;
82  long x, y;
83  amg_point *pointx, *pointy;
84  unsigned int c_points = Pointvector[level].get_cpoints();
85 
86  // Setup Prolongation/Interpolation matrix
87  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()),c_points);
88  P[level].clear();
89 
90  // Assign indices to C points
91  Pointvector[level].build_index();
92 
93  // Direct Interpolation (Yang, p.14)
94 #ifdef VIENNACL_WITH_OPENMP
95  #pragma omp parallel for private (pointx,pointy,row_sum,c_sum,temp_res,y,x,diag)
96 #endif
97  for (x=0; x < static_cast<long>(Pointvector[level].size()); ++x)
98  {
99  pointx = Pointvector[level][x];
100  /*if (A[level](x,x) > 0)
101  diag_sign = 1;
102  else
103  diag_sign = -1;*/
104 
105  // When the current line corresponds to a C point then the diagonal coefficient is 1 and the rest 0
106  if (pointx->is_cpoint())
107  P[level](x,pointx->get_coarse_index()) = 1;
108 
109  // When the current line corresponds to a F point then the diagonal is 0 and the rest has to be computed (Yang, p.14)
110  if (pointx->is_fpoint())
111  {
112  // Jump to row x
113  InternalRowIterator row_iter = A[level].begin1();
114  row_iter += x;
115 
116  // Row sum of coefficients (without diagonal) and sum of influencing C point coefficients has to be computed
117  row_sum = c_sum = diag = 0;
118  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
119  {
120  y = static_cast<long>(col_iter.index2());
121  if (x == y)// || *col_iter * diag_sign > 0)
122  {
123  diag += *col_iter;
124  continue;
125  }
126 
127  // Sum all other coefficients in line x
128  row_sum += *col_iter;
129 
130  pointy = Pointvector[level][y];
131  // Sum all coefficients that correspond to a strongly influencing C point
132  if (pointy->is_cpoint())
133  if (pointx->is_influencing(pointy))
134  c_sum += *col_iter;
135  }
136  temp_res = -row_sum/(c_sum*diag);
137 
138  // Iterate over all strongly influencing points of point x
139  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
140  {
141  pointy = *iter;
142  // The value is only non-zero for columns that correspond to a C point
143  if (pointy->is_cpoint())
144  {
145  if (temp_res != 0)
146  P[level](x, pointy->get_coarse_index()) = temp_res * A[level](x,pointy->get_index());
147  }
148  }
149 
150  //Truncate interpolation if chosen
151  if (tag.get_interpolweight() != 0)
152  amg_truncate_row(P[level], x, tag);
153  }
154  }
155 
156  // P test
157  //test_interpolation(A[level], P[level], Pointvector[level]);
158 
159  #ifdef VIENNACL_AMG_DEBUG
160  std::cout << "Prolongation Matrix:" << std::endl;
161  printmatrix (P[level]);
162  #endif
163  }
164 
172  template <typename InternalType1, typename InternalType2>
173  void amg_interpol_classic(unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
174  {
175  typedef typename InternalType1::value_type SparseMatrixType;
176  //typedef typename InternalType2::value_type PointVectorType;
177  typedef typename SparseMatrixType::value_type ScalarType;
178  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
179  typedef typename SparseMatrixType::iterator2 InternalColIterator;
180 
181  ScalarType temp_res;
182  ScalarType weak_sum, strong_sum;
183  int diag_sign;
185  amg_point *pointx, *pointy, *pointk, *pointm;
186  long x, y, k, m;
187 
188  unsigned int c_points = Pointvector[level].get_cpoints();
189 
190  // Setup Prolongation/Interpolation matrix
191  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
192  P[level].clear();
193 
194  // Assign indices to C points
195  Pointvector[level].build_index();
196 
197  // Classical Interpolation (Yang, p.13-14)
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)
200 #endif
201  for (x=0; x < static_cast<long>(Pointvector[level].size()); ++x)
202  {
203  pointx = Pointvector[level][x];
204  if (A[level](x,x) > 0)
205  diag_sign = 1;
206  else
207  diag_sign = -1;
208 
209  // When the current line corresponds to a C point then the diagonal coefficient is 1 and the rest 0
210  if (pointx->is_cpoint())
211  P[level](x,pointx->get_coarse_index()) = 1;
212 
213  // When the current line corresponds to a F point then the diagonal is 0 and the rest has to be computed (Yang, p.14)
214  if (pointx->is_fpoint())
215  {
216  // Jump to row x
217  InternalRowIterator row_iter = A[level].begin1();
218  row_iter += x;
219 
220  weak_sum = 0;
221  c_sum_row = amg_sparsevector<ScalarType>(static_cast<unsigned int>(A[level].size1()));
222  c_sum_row.clear();
223  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
224  {
225  k = static_cast<unsigned int>(col_iter.index2());
226  pointk = Pointvector[level][k];
227 
228  // Sum of weakly influencing neighbors + diagonal coefficient
229  if (x == k || !pointx->is_influencing(pointk))// || *col_iter * diag_sign > 0)
230  {
231  weak_sum += *col_iter;
232  continue;
233  }
234 
235  // Sums of coefficients in row k (strongly influening F neighbors) of C point neighbors of x are calculated
236  if (pointk->is_fpoint() && pointx->is_influencing(pointk))
237  {
238  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
239  {
240  pointm = *iter;
241  m = pointm->get_index();
242 
243  if (pointm->is_cpoint())
244  // Only use coefficients that have opposite sign of diagonal.
245  if (A[level](k,m) * diag_sign < 0)
246  c_sum_row[k] += A[level](k,m);
247  }
248  continue;
249  }
250  }
251 
252  // Iterate over all strongly influencing points of point x
253  for (amg_point::iterator iter = pointx->begin_influencing(); iter != pointx->end_influencing(); ++iter)
254  {
255  pointy = *iter;
256  y = pointy->get_index();
257 
258  // The value is only non-zero for columns that correspond to a C point
259  if (pointy->is_cpoint())
260  {
261  strong_sum = 0;
262  // Calculate term for strongly influencing F neighbors
263  for (typename amg_sparsevector<ScalarType>::iterator iter2 = c_sum_row.begin(); iter2 != c_sum_row.end(); ++iter2)
264  {
265  k = iter2.index();
266  // Only use coefficients that have opposite sign of diagonal.
267  if (A[level](k,y) * diag_sign < 0)
268  strong_sum += (A[level](x,k) * A[level](k,y)) / (*iter2);
269  }
270 
271  // Calculate coefficient
272  temp_res = - (A[level](x,y) + strong_sum) / (weak_sum);
273  if (temp_res != 0)
274  P[level](x,pointy->get_coarse_index()) = temp_res;
275  }
276  }
277 
278  //Truncate iteration if chosen
279  if (tag.get_interpolweight() != 0)
280  amg_truncate_row(P[level], x, tag);
281  }
282  }
283 
284  #ifdef VIENNACL_AMG_DEBUG
285  std::cout << "Prolongation Matrix:" << std::endl;
286  printmatrix (P[level]);
287  #endif
288  }
289 
296  template <typename SparseMatrixType>
297  void amg_truncate_row(SparseMatrixType & P, unsigned int row, amg_tag & tag)
298  {
299  typedef typename SparseMatrixType::value_type ScalarType;
300  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
301  typedef typename SparseMatrixType::iterator2 InternalColIterator;
302 
303  ScalarType row_max, row_min, row_sum_pos, row_sum_neg, row_sum_pos_scale, row_sum_neg_scale;
304 
305  InternalRowIterator row_iter = P.begin1();
306  row_iter += row;
307 
308  row_max = 0;
309  row_min = 0;
310  row_sum_pos = 0;
311  row_sum_neg = 0;
312 
313  // Truncate interpolation by making values to zero that are a lot smaller than the biggest value in a row
314  // Determine max entry and sum of row (seperately for negative and positive entries)
315  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
316  {
317  if (*col_iter > row_max)
318  row_max = *col_iter;
319  if (*col_iter < row_min)
320  row_min = *col_iter;
321  if (*col_iter > 0)
322  row_sum_pos += *col_iter;
323  if (*col_iter < 0)
324  row_sum_neg += *col_iter;
325  }
326 
327  row_sum_pos_scale = row_sum_pos;
328  row_sum_neg_scale = row_sum_neg;
329 
330  // Make certain values to zero (seperately for negative and positive entries)
331  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
332  {
333  if (*col_iter > 0 && *col_iter < tag.get_interpolweight() * row_max)
334  {
335  row_sum_pos_scale -= *col_iter;
336  *col_iter = 0;
337  }
338  if (*col_iter < 0 && *col_iter > tag.get_interpolweight() * row_min)
339  {
340  row_sum_pos_scale -= *col_iter;
341  *col_iter = 0;
342  }
343  }
344 
345  // Scale remaining values such that row sum is unchanged
346  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
347  {
348  if (*col_iter > 0)
349  *col_iter = *col_iter *(row_sum_pos/row_sum_pos_scale);
350  if (*col_iter < 0)
351  *col_iter = *col_iter *(row_sum_neg/row_sum_neg_scale);
352  }
353  }
354 
361  template <typename InternalType1, typename InternalType2>
362  void amg_interpol_ag(unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag)
363  {
364  typedef typename InternalType1::value_type SparseMatrixType;
365  //typedef typename InternalType2::value_type PointVectorType;
366  //typedef typename SparseMatrixType::value_type ScalarType;
367  //typedef typename SparseMatrixType::iterator1 InternalRowIterator;
368  //typedef typename SparseMatrixType::iterator2 InternalColIterator;
369 
370  long x;
371  amg_point *pointx, *pointy;
372  unsigned int c_points = Pointvector[level].get_cpoints();
373 
374  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
375  P[level].clear();
376 
377  // Assign indices to C points
378  Pointvector[level].build_index();
379 
380  // Set prolongation such that F point is interpolated (weight=1) by the aggregate it belongs to (Vanek et al p.6)
381 #ifdef VIENNACL_WITH_OPENMP
382  #pragma omp parallel for private (x,pointx)
383 #endif
384  for (x=0; x<static_cast<long>(Pointvector[level].size()); ++x)
385  {
386  pointx = Pointvector[level][x];
387  pointy = Pointvector[level][pointx->get_aggregate()];
388  // Point x belongs to aggregate y.
389  P[level](x,pointy->get_coarse_index()) = 1;
390  }
391 
392  #ifdef VIENNACL_AMG_DEBUG
393  std::cout << "Aggregation based Prolongation:" << std::endl;
394  printmatrix(P[level]);
395  #endif
396  }
397 
405  template <typename InternalType1, typename InternalType2>
406  void amg_interpol_sa(unsigned int level, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
407  {
408  typedef typename InternalType1::value_type SparseMatrixType;
409  //typedef typename InternalType2::value_type PointVectorType;
410  typedef typename SparseMatrixType::value_type ScalarType;
411  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
412  typedef typename SparseMatrixType::iterator2 InternalColIterator;
413 
414  long x,y;
415  ScalarType diag = 0;
416  unsigned int c_points = Pointvector[level].get_cpoints();
417 
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()));
420  Jacobi.clear();
421  P[level] = SparseMatrixType(static_cast<unsigned int>(A[level].size1()), c_points);
422  P[level].clear();
423 
424  // Build Jacobi Matrix via filtered A matrix (Vanek et al. p.6)
425 #ifdef VIENNACL_WITH_OPENMP
426  #pragma omp parallel for private (x,y,diag)
427 #endif
428  for (x=0; x<static_cast<long>(A[level].size1()); ++x)
429  {
430  diag = 0;
431  InternalRowIterator row_iter = A[level].begin1();
432  row_iter += x;
433  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
434  {
435  y = static_cast<long>(col_iter.index2());
436  // Determine the structure of the Jacobi matrix by using a filtered matrix of A:
437  // The diagonal consists of the diagonal coefficient minus all coefficients of points not in the neighborhood of x.
438  // All other coefficients are the same as in A.
439  // Already use Jacobi matrix to save filtered A matrix to speed up computation.
440  if (x == y)
441  diag += *col_iter;
442  else if (!Pointvector[level][x]->is_influencing(Pointvector[level][y]))
443  diag += -*col_iter;
444  else
445  Jacobi (x,y) = *col_iter;
446  }
447  InternalRowIterator row_iter2 = Jacobi.begin1();
448  row_iter2 += x;
449  // Traverse through filtered A matrix and compute the Jacobi filtering
450  for (InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
451  {
452  *col_iter2 = - static_cast<ScalarType>(tag.get_interpolweight())/diag * *col_iter2;
453  }
454  // Diagonal can be computed seperately.
455  Jacobi (x,x) = 1 - static_cast<ScalarType>(tag.get_interpolweight());
456  }
457 
458  #ifdef VIENNACL_AMG_DEBUG
459  std::cout << "Jacobi Matrix:" << std::endl;
460  printmatrix(Jacobi);
461  #endif
462 
463  // Use AG interpolation as tentative prolongation
464  amg_interpol_ag(level, A, P_tentative, Pointvector, tag);
465 
466  #ifdef VIENNACL_AMG_DEBUG
467  std::cout << "Tentative Prolongation:" << std::endl;
468  printmatrix(P_tentative[level]);
469  #endif
470 
471  // Multiply Jacobi matrix with tentative prolongation to get actual prolongation
472  amg_mat_prod(Jacobi,P_tentative[level],P[level]);
473 
474  #ifdef VIENNACL_AMG_DEBUG
475  std::cout << "Prolongation Matrix:" << std::endl;
476  printmatrix (P[level]);
477  #endif
478  }
479  } //namespace amg
480  }
481  }
482 }
483 
484 #endif
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
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