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