ViennaCL - The Vienna Computing Library  1.5.1
amg.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_AMG_HPP_
2 #define VIENNACL_LINALG_AMG_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 
27 #include <boost/numeric/ublas/matrix.hpp>
28 #include <boost/numeric/ublas/lu.hpp>
29 #include <boost/numeric/ublas/operation.hpp>
30 #include <boost/numeric/ublas/vector_proxy.hpp>
31 #include <boost/numeric/ublas/matrix_proxy.hpp>
32 #include <boost/numeric/ublas/vector.hpp>
33 #include <boost/numeric/ublas/triangular.hpp>
34 #include <vector>
35 #include <cmath>
36 #include "viennacl/forwards.h"
37 #include "viennacl/tools/tools.hpp"
38 #include "viennacl/linalg/prod.hpp"
40 
44 
45 #include <map>
46 
47 #ifdef VIENNACL_WITH_OPENMP
48  #include <omp.h>
49 #endif
50 
52 
53 #define VIENNACL_AMG_COARSE_LIMIT 50
54 #define VIENNACL_AMG_MAX_LEVELS 100
55 
56 namespace viennacl
57 {
58  namespace linalg
59  {
61 
62 
63 
71  template <typename InternalType1, typename InternalType2>
72  void amg_setup(InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
73  {
74  typedef typename InternalType2::value_type PointVectorType;
75 
76  unsigned int i, iterations, c_points, f_points;
78 
79  // Set number of iterations. If automatic coarse grid construction is chosen (0), then set a maximum size and stop during the process.
80  iterations = tag.get_coarselevels();
81  if (iterations == 0)
82  iterations = VIENNACL_AMG_MAX_LEVELS;
83 
84  // For parallel coarsenings build data structures (number of threads set automatically).
86  Slicing.init(iterations);
87 
88  for (i=0; i<iterations; ++i)
89  {
90  // Initialize Pointvector on level i and construct points.
91  Pointvector[i] = PointVectorType(static_cast<unsigned int>(A[i].size1()));
92  Pointvector[i].init_points();
93 
94  // Construct C and F points on coarse level (i is fine level, i+1 coarse level).
95  detail::amg::amg_coarse (i, A, Pointvector, Slicing, tag);
96 
97  // Calculate number of C and F points on level i.
98  c_points = Pointvector[i].get_cpoints();
99  f_points = Pointvector[i].get_fpoints();
100 
101  #if defined (VIENNACL_AMG_DEBUG) //or defined(VIENNACL_AMG_DEBUGBENCH)
102  std::cout << "Level " << i << ": ";
103  std::cout << "No of C points = " << c_points << ", ";
104  std::cout << "No of F points = " << f_points << std::endl;
105  #endif
106 
107  // Stop routine when the maximal coarse level is found (no C or F point). Coarsest level is level i.
108  if (c_points == 0 || f_points == 0)
109  break;
110 
111  // Construct interpolation matrix for level i.
112  detail::amg::amg_interpol (i, A, P, Pointvector, tag);
113 
114  // Compute coarse grid operator (A[i+1] = R * A[i] * P) with R = trans(P).
115  detail::amg::amg_galerkin_prod(A[i], P[i], A[i+1]);
116 
117  // Test triple matrix product. Very slow for large matrix sizes (ublas).
118  // test_triplematprod(A[i],P[i],A[i+1]);
119 
120  Pointvector[i].delete_points();
121 
122  #ifdef VIENNACL_AMG_DEBUG
123  std::cout << "Coarse Grid Operator Matrix:" << std::endl;
124  printmatrix (A[i+1]);
125  #endif
126 
127  // If Limit of coarse points is reached then stop. Coarsest level is level i+1.
128  if (tag.get_coarselevels() == 0 && c_points <= VIENNACL_AMG_COARSE_LIMIT)
129  {
130  tag.set_coarselevels(i+1);
131  return;
132  }
133  }
134  tag.set_coarselevels(i);
135  }
136 
145  template <typename MatrixType, typename InternalType1, typename InternalType2>
146  void amg_init(MatrixType const & mat, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
147  {
148  //typedef typename MatrixType::value_type ScalarType;
149  typedef typename InternalType1::value_type SparseMatrixType;
150 
151  if (tag.get_coarselevels() > 0)
152  {
153  A.resize(tag.get_coarselevels()+1);
154  P.resize(tag.get_coarselevels());
155  Pointvector.resize(tag.get_coarselevels());
156  }
157  else
158  {
159  A.resize(VIENNACL_AMG_MAX_LEVELS+1);
160  P.resize(VIENNACL_AMG_MAX_LEVELS);
161  Pointvector.resize(VIENNACL_AMG_MAX_LEVELS);
162  }
163 
164  // Insert operator matrix as operator for finest level.
165  SparseMatrixType A0 (mat);
166  A.insert_element (0, A0);
167  }
168 
178  template <typename InternalType1, typename InternalType2>
179  void amg_transform_cpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup, amg_tag & tag)
180  {
181  //typedef typename InternalType1::value_type MatrixType;
182 
183  // Resize internal data structures to actual size.
184  A.resize(tag.get_coarselevels()+1);
185  P.resize(tag.get_coarselevels());
186  R.resize(tag.get_coarselevels());
187 
188  // Transform into matrix type.
189  for (unsigned int i=0; i<tag.get_coarselevels()+1; ++i)
190  {
191  A[i].resize(A_setup[i].size1(),A_setup[i].size2(),false);
192  A[i] = A_setup[i];
193  }
194  for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
195  {
196  P[i].resize(P_setup[i].size1(),P_setup[i].size2(),false);
197  P[i] = P_setup[i];
198  }
199  for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
200  {
201  R[i].resize(P_setup[i].size2(),P_setup[i].size1(),false);
202  P_setup[i].set_trans(true);
203  R[i] = P_setup[i];
204  P_setup[i].set_trans(false);
205  }
206  }
207 
218  template <typename InternalType1, typename InternalType2>
219  void amg_transform_gpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup, amg_tag & tag, viennacl::context ctx)
220  {
221  // Resize internal data structures to actual size.
222  A.resize(tag.get_coarselevels()+1);
223  P.resize(tag.get_coarselevels());
224  R.resize(tag.get_coarselevels());
225 
226  // Copy to GPU using the internal sparse matrix structure: std::vector<std::map>.
227  for (unsigned int i=0; i<tag.get_coarselevels()+1; ++i)
228  {
230  //A[i].resize(A_setup[i].size1(),A_setup[i].size2(),false);
231  viennacl::copy(*(A_setup[i].get_internal_pointer()),A[i]);
232  }
233  for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
234  {
236  //P[i].resize(P_setup[i].size1(),P_setup[i].size2(),false);
237  viennacl::copy(*(P_setup[i].get_internal_pointer()),P[i]);
238  //viennacl::copy((boost::numeric::ublas::compressed_matrix<ScalarType>)P_setup[i],P[i]);
239  }
240  for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
241  {
243  //R[i].resize(P_setup[i].size2(),P_setup[i].size1(),false);
244  P_setup[i].set_trans(true);
245  viennacl::copy(*(P_setup[i].get_internal_pointer()),R[i]);
246  P_setup[i].set_trans(false);
247  }
248  }
249 
258  template <typename InternalVectorType, typename SparseMatrixType>
259  void amg_setup_apply (InternalVectorType & result, InternalVectorType & rhs, InternalVectorType & residual, SparseMatrixType const & A, amg_tag const & tag)
260  {
261  typedef typename InternalVectorType::value_type VectorType;
262 
263  result.resize(tag.get_coarselevels()+1);
264  rhs.resize(tag.get_coarselevels()+1);
265  residual.resize(tag.get_coarselevels());
266 
267  for (unsigned int level=0; level < tag.get_coarselevels()+1; ++level)
268  {
269  result[level] = VectorType(A[level].size1());
270  result[level].clear();
271  rhs[level] = VectorType(A[level].size1());
272  rhs[level].clear();
273  }
274  for (unsigned int level=0; level < tag.get_coarselevels(); ++level)
275  {
276  residual[level] = VectorType(A[level].size1());
277  residual[level].clear();
278  }
279  }
280 
281 
291  template <typename InternalVectorType, typename SparseMatrixType>
292  void amg_setup_apply (InternalVectorType & result, InternalVectorType & rhs, InternalVectorType & residual, SparseMatrixType const & A, amg_tag const & tag, viennacl::context ctx)
293  {
294  typedef typename InternalVectorType::value_type VectorType;
295 
296  result.resize(tag.get_coarselevels()+1);
297  rhs.resize(tag.get_coarselevels()+1);
298  residual.resize(tag.get_coarselevels());
299 
300  for (unsigned int level=0; level < tag.get_coarselevels()+1; ++level)
301  {
302  result[level] = VectorType(A[level].size1(), ctx);
303  rhs[level] = VectorType(A[level].size1(), ctx);
304  }
305  for (unsigned int level=0; level < tag.get_coarselevels(); ++level)
306  {
307  residual[level] = VectorType(A[level].size1(), ctx);
308  }
309  }
310 
311 
319  template <typename ScalarType, typename SparseMatrixType>
320  void amg_lu(boost::numeric::ublas::compressed_matrix<ScalarType> & op, boost::numeric::ublas::permutation_matrix<> & Permutation, SparseMatrixType const & A)
321  {
322  typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
323  typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
324 
325  // Copy to operator matrix. Needed
326  op.resize(A.size1(),A.size2(),false);
327  for (ConstRowIterator row_iter = A.begin1(); row_iter != A.end1(); ++row_iter)
328  for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
329  op (col_iter.index1(), col_iter.index2()) = *col_iter;
330 
331  // Permutation matrix has to be reinitialized with actual size. Do not clear() or resize()!
332  Permutation = boost::numeric::ublas::permutation_matrix<> (op.size1());
333  boost::numeric::ublas::lu_factorize(op,Permutation);
334  }
335 
338  template <typename MatrixType>
340  {
341  typedef typename MatrixType::value_type ScalarType;
342  typedef boost::numeric::ublas::vector<ScalarType> VectorType;
345 
350 
351  boost::numeric::ublas::vector <SparseMatrixType> A_setup;
352  boost::numeric::ublas::vector <SparseMatrixType> P_setup;
353  boost::numeric::ublas::vector <MatrixType> A;
354  boost::numeric::ublas::vector <MatrixType> P;
355  boost::numeric::ublas::vector <MatrixType> R;
356  boost::numeric::ublas::vector <PointVectorType> Pointvector;
357 
358  mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
359  mutable boost::numeric::ublas::permutation_matrix<> Permutation;
360 
361  mutable boost::numeric::ublas::vector <VectorType> result;
362  mutable boost::numeric::ublas::vector <VectorType> rhs;
363  mutable boost::numeric::ublas::vector <VectorType> residual;
364 
365  mutable bool done_init_apply;
366 
367  amg_tag tag_;
368  public:
369 
370  amg_precond(): Permutation(0) {}
376  amg_precond(MatrixType const & mat, amg_tag const & tag): Permutation(0)
377  {
378  tag_ = tag;
379  // Initialize data structures.
380  amg_init (mat,A_setup,P_setup,Pointvector,tag_);
381 
382  done_init_apply = false;
383  }
384 
387  void setup()
388  {
389  // Start setup phase.
390  amg_setup(A_setup,P_setup,Pointvector,tag_);
391  // Transform to CPU-Matrixtype for precondition phase.
392  amg_transform_cpu(A,P,R,A_setup,P_setup,tag_);
393 
394  done_init_apply = false;
395  }
396 
401  void init_apply() const
402  {
403  // Setup precondition phase (Data structures).
404  amg_setup_apply(result,rhs,residual,A_setup,tag_);
405  // Do LU factorization for direct solve.
406  amg_lu(op,Permutation,A_setup[tag_.get_coarselevels()]);
407 
408  done_init_apply = true;
409  }
410 
416  template <typename VectorType>
417  ScalarType calc_complexity(VectorType & avgstencil)
418  {
419  avgstencil = VectorType (tag_.get_coarselevels()+1);
420  unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
421 
422  for (unsigned int level=0; level < tag_.get_coarselevels()+1; ++level)
423  {
424  level_coefficients = 0;
425  for (InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
426  {
427  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
428  {
429  if (level == 0)
430  systemmat_nonzero++;
431  nonzero++;
432  level_coefficients++;
433  }
434  }
435  avgstencil[level] = level_coefficients/static_cast<ScalarType>(A_setup[level].size1());
436  }
437  return nonzero/static_cast<ScalarType>(systemmat_nonzero);
438  }
439 
444  template <typename VectorType>
445  void apply(VectorType & vec) const
446  {
447  // Build data structures and do lu factorization before first iteration step.
448  if (!done_init_apply)
449  init_apply();
450 
451  int level;
452 
453  // Precondition operation (Yang, p.3)
454  rhs[0] = vec;
455  for (level=0; level <static_cast<int>(tag_.get_coarselevels()); level++)
456  {
457  result[level].clear();
458 
459  // Apply Smoother presmooth_ times.
460  smooth_jacobi (level, tag_.get_presmooth(), result[level], rhs[level]);
461 
462  #ifdef VIENNACL_AMG_DEBUG
463  std::cout << "After presmooth:" << std::endl;
464  printvector(result[level]);
465  #endif
466 
467  // Compute residual.
468  residual[level] = rhs[level] - boost::numeric::ublas::prod (A[level],result[level]);
469 
470  #ifdef VIENNACL_AMG_DEBUG
471  std::cout << "Residual:" << std::endl;
472  printvector(residual[level]);
473  #endif
474 
475  // Restrict to coarse level. Restricted residual is RHS of coarse level.
476  rhs[level+1] = boost::numeric::ublas::prod (R[level],residual[level]);
477 
478  #ifdef VIENNACL_AMG_DEBUG
479  std::cout << "Restricted Residual: " << std::endl;
480  printvector(rhs[level+1]);
481  #endif
482  }
483 
484  // On highest level use direct solve to solve equation.
485  result[level] = rhs[level];
486  boost::numeric::ublas::lu_substitute(op,Permutation,result[level]);
487 
488  #ifdef VIENNACL_AMG_DEBUG
489  std::cout << "After direct solve: " << std::endl;
490  printvector (result[level]);
491  #endif
492 
493  for (level=tag_.get_coarselevels()-1; level >= 0; level--)
494  {
495  #ifdef VIENNACL_AMG_DEBUG
496  std::cout << "Coarse Error: " << std::endl;
497  printvector(result[level+1]);
498  #endif
499 
500  // Interpolate error to fine level. Correct solution by adding error.
501  result[level] += boost::numeric::ublas::prod (P[level], result[level+1]);
502 
503  #ifdef VIENNACL_AMG_DEBUG
504  std::cout << "Corrected Result: " << std::endl;
505  printvector (result[level]);
506  #endif
507 
508  // Apply Smoother postsmooth_ times.
509  smooth_jacobi (level, tag_.get_postsmooth(), result[level], rhs[level]);
510 
511  #ifdef VIENNACL_AMG_DEBUG
512  std::cout << "After postsmooth: " << std::endl;
513  printvector (result[level]);
514  #endif
515  }
516  vec = result[0];
517  }
518 
525  template <typename VectorType>
526  void smooth_jacobi(int level, int const iterations, VectorType & x, VectorType const & rhs) const
527  {
528  VectorType old_result (x.size());
529  long index;
530  ScalarType sum = 0, diag = 1;
531 
532  for (int i=0; i<iterations; ++i)
533  {
534  old_result = x;
535  x.clear();
536 #ifdef VIENNACL_WITH_OPENMP
537  #pragma omp parallel for private (sum,diag) shared (rhs,x)
538 #endif
539  for (index=0; index < static_cast<long>(A_setup[level].size1()); ++index)
540  {
541  InternalConstRowIterator row_iter = A_setup[level].begin1();
542  row_iter += index;
543  sum = 0;
544  diag = 1;
545  for (InternalConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
546  {
547  if (col_iter.index1() == col_iter.index2())
548  diag = *col_iter;
549  else
550  sum += *col_iter * old_result[col_iter.index2()];
551  }
552  x[index]= static_cast<ScalarType>(tag_.get_jacobiweight()) * (rhs[index] - sum) / diag + (1-static_cast<ScalarType>(tag_.get_jacobiweight())) * old_result[index];
553  }
554  }
555  }
556 
557  amg_tag & tag() { return tag_; }
558  };
559 
564  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
565  class amg_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
566  {
571 
576 
577  boost::numeric::ublas::vector <SparseMatrixType> A_setup;
578  boost::numeric::ublas::vector <SparseMatrixType> P_setup;
579  boost::numeric::ublas::vector <MatrixType> A;
580  boost::numeric::ublas::vector <MatrixType> P;
581  boost::numeric::ublas::vector <MatrixType> R;
582  boost::numeric::ublas::vector <PointVectorType> Pointvector;
583 
584  mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
585  mutable boost::numeric::ublas::permutation_matrix<> Permutation;
586 
587  mutable boost::numeric::ublas::vector <VectorType> result;
588  mutable boost::numeric::ublas::vector <VectorType> rhs;
589  mutable boost::numeric::ublas::vector <VectorType> residual;
590 
591  viennacl::context ctx_;
592 
593  mutable bool done_init_apply;
594 
595  amg_tag tag_;
596 
597  public:
598 
599  amg_precond(): Permutation(0) {}
600 
606  amg_precond(compressed_matrix<ScalarType, MAT_ALIGNMENT> const & mat, amg_tag const & tag): Permutation(0), ctx_(viennacl::traits::context(mat))
607  {
608  tag_ = tag;
609 
610  // Copy to CPU. Internal structure of sparse matrix is used for copy operation.
611  std::vector<std::map<unsigned int, ScalarType> > mat2 = std::vector<std::map<unsigned int, ScalarType> >(mat.size1());
612  viennacl::copy(mat, mat2);
613 
614  // Initialize data structures.
615  amg_init (mat2,A_setup,P_setup,Pointvector,tag_);
616 
617  done_init_apply = false;
618  }
619 
622  void setup()
623  {
624  // Start setup phase.
625  amg_setup(A_setup,P_setup,Pointvector, tag_);
626  // Transform to GPU-Matrixtype for precondition phase.
627  amg_transform_gpu(A,P,R,A_setup,P_setup, tag_, ctx_);
628 
629  done_init_apply = false;
630  }
631 
636  void init_apply() const
637  {
638  // Setup precondition phase (Data structures).
639  amg_setup_apply(result,rhs,residual,A_setup,tag_, ctx_);
640  // Do LU factorization for direct solve.
641  amg_lu(op,Permutation,A_setup[tag_.get_coarselevels()]);
642 
643  done_init_apply = true;
644  }
645 
651  template <typename VectorType>
652  ScalarType calc_complexity(VectorType & avgstencil)
653  {
654  avgstencil = VectorType (tag_.get_coarselevels()+1);
655  unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
656 
657  for (unsigned int level=0; level < tag_.get_coarselevels()+1; ++level)
658  {
659  level_coefficients = 0;
660  for (InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
661  {
662  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
663  {
664  if (level == 0)
665  systemmat_nonzero++;
666  nonzero++;
667  level_coefficients++;
668  }
669  }
670  avgstencil[level] = level_coefficients/(double)A[level].size1();
671  }
672  return nonzero/static_cast<double>(systemmat_nonzero);
673  }
674 
679  template <typename VectorType>
680  void apply(VectorType & vec) const
681  {
682  if (!done_init_apply)
683  init_apply();
684 
685  int level;
686 
687  // Precondition operation (Yang, p.3).
688  rhs[0] = vec;
689  for (level=0; level <static_cast<int>(tag_.get_coarselevels()); level++)
690  {
691  result[level].clear();
692 
693  // Apply Smoother presmooth_ times.
694  smooth_jacobi (level, tag_.get_presmooth(), result[level], rhs[level]);
695 
696  #ifdef VIENNACL_AMG_DEBUG
697  std::cout << "After presmooth: " << std::endl;
698  printvector(result[level]);
699  #endif
700 
701  // Compute residual.
702  //residual[level] = rhs[level] - viennacl::linalg::prod (A[level],result[level]);
703  residual[level] = viennacl::linalg::prod (A[level],result[level]);
704  residual[level] = rhs[level] - residual[level];
705 
706  #ifdef VIENNACL_AMG_DEBUG
707  std::cout << "Residual: " << std::endl;
708  printvector(residual[level]);
709  #endif
710 
711  // Restrict to coarse level. Result is RHS of coarse level equation.
712  //residual_coarse[level] = viennacl::linalg::prod(R[level],residual[level]);
713  rhs[level+1] = viennacl::linalg::prod(R[level],residual[level]);
714 
715  #ifdef VIENNACL_AMG_DEBUG
716  std::cout << "Restricted Residual: " << std::endl;
717  printvector(rhs[level+1]);
718  #endif
719  }
720 
721  // On highest level use direct solve to solve equation (on the CPU)
722  //TODO: Use GPU direct solve!
723  result[level] = rhs[level];
724  boost::numeric::ublas::vector <ScalarType> result_cpu (result[level].size());
725 
726  copy (result[level],result_cpu);
727  boost::numeric::ublas::lu_substitute(op,Permutation,result_cpu);
728  copy (result_cpu, result[level]);
729 
730  #ifdef VIENNACL_AMG_DEBUG
731  std::cout << "After direct solve: " << std::endl;
732  printvector (result[level]);
733  #endif
734 
735  for (level=tag_.get_coarselevels()-1; level >= 0; level--)
736  {
737  #ifdef VIENNACL_AMG_DEBUG
738  std::cout << "Coarse Error: " << std::endl;
739  printvector(result[level+1]);
740  #endif
741 
742  // Interpolate error to fine level and correct solution.
743  result[level] += viennacl::linalg::prod(P[level],result[level+1]);
744 
745  #ifdef VIENNACL_AMG_DEBUG
746  std::cout << "Corrected Result: " << std::endl;
747  printvector (result[level]);
748  #endif
749 
750  // Apply Smoother postsmooth_ times.
751  smooth_jacobi (level, tag_.get_postsmooth(), result[level], rhs[level]);
752 
753  #ifdef VIENNACL_AMG_DEBUG
754  std::cout << "After postsmooth: " << std::endl;
755  printvector (result[level]);
756  #endif
757  }
758  vec = result[0];
759  }
760 
767  template <typename VectorType>
768  void smooth_jacobi(int level, unsigned int iterations, VectorType & x, VectorType const & rhs) const
769  {
770  VectorType old_result = x;
771 
772  viennacl::ocl::context & ctx = const_cast<viennacl::ocl::context &>(viennacl::traits::opencl_handle(x).context());
775 
776  for (unsigned int i=0; i<iterations; ++i)
777  {
778  if (i > 0)
779  old_result = x;
780  x.clear();
781  viennacl::ocl::enqueue(k(A[level].handle1().opencl_handle(), A[level].handle2().opencl_handle(), A[level].handle().opencl_handle(),
782  static_cast<ScalarType>(tag_.get_jacobiweight()),
783  viennacl::traits::opencl_handle(old_result),
784  viennacl::traits::opencl_handle(x),
785  viennacl::traits::opencl_handle(rhs),
786  static_cast<cl_uint>(rhs.size())));
787 
788  }
789  }
790 
791  amg_tag & tag() { return tag_; }
792  };
793 
794  }
795 }
796 
797 
798 
799 #endif
800 
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
viennacl::ocl::kernel & get_kernel(std::string const &program_name, std::string const &kernel_name)
Convenience function for retrieving the kernel of a program directly from the context.
Definition: context.hpp:470
Debug functionality for AMG. To be removed.
double get_jacobiweight() const
Definition: amg_base.hpp:102
void set_coarselevels(int coarselevels)
Definition: amg_base.hpp:110
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:859
void amg_lu(boost::numeric::ublas::compressed_matrix< ScalarType > &op, boost::numeric::ublas::permutation_matrix<> &Permutation, SparseMatrixType const &A)
Pre-compute LU factorization for direct solve (ublas library).
Definition: amg.hpp:320
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Represents an OpenCL kernel within ViennaCL.
Definition: kernel.hpp:59
amg_precond(compressed_matrix< ScalarType, MAT_ALIGNMENT > const &mat, amg_tag const &tag)
The constructor. Builds data structures.
Definition: amg.hpp:606
Various little tools used here and there in ViennaCL.
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
#define VIENNACL_AMG_COARSE_RS3
Definition: amg_base.hpp:44
Manages an OpenCL context and provides the respective convenience functions for creating buffers...
Definition: context.hpp:51
void prod(const T1 &A, bool transposed_A, const T2 &B, bool transposed_B, T3 &C, ScalarType alpha, ScalarType beta)
Definition: matrix_operations.hpp:2305
void init(unsigned int levels, unsigned int threads=0)
Definition: amg_base.hpp:1183
unsigned int get_coarse() const
Definition: amg_base.hpp:90
void amg_init(MatrixType const &mat, InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Initialize AMG preconditioner.
Definition: amg.hpp:146
This file provides the forward declarations for the main types used within ViennaCL.
sparse_matrix_adapted_iterator< SCALARTYPE, SizeType,!is_iterator1 > begin() const
Definition: adapter.hpp:319
const_sparse_matrix_adapted_iterator< SCALARTYPE, SizeType,!is_iterator1, true > end() const
Definition: adapter.hpp:156
unsigned int get_postsmooth() const
Definition: amg_base.hpp:108
void printmatrix(MatrixType &, int)
Definition: amg_debug.hpp:77
void init_apply() const
Prepare data structures for preconditioning: Build data structures for precondition phase...
Definition: amg.hpp:401
void amg_transform_cpu(InternalType1 &A, InternalType1 &P, InternalType1 &R, InternalType2 &A_setup, InternalType2 &P_setup, amg_tag &tag)
Save operators after setup phase for CPU computation.
Definition: amg.hpp:179
#define VIENNACL_AMG_COARSE_LIMIT
Definition: amg.hpp:53
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
void apply(VectorType &vec) const
Precondition Operation.
Definition: amg.hpp:680
#define VIENNACL_AMG_COARSE_RS0
Definition: amg_base.hpp:43
ScalarType calc_complexity(VectorType &avgstencil)
Returns complexity measures.
Definition: amg.hpp:417
void amg_setup(InternalType1 &A, InternalType1 &P, InternalType2 &Pointvector, amg_tag &tag)
Setup AMG preconditioner.
Definition: amg.hpp:72
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
A const iterator for sparse matrices of type std::vector<std::map<SizeType, SCALARTYPE> > ...
Definition: adapter.hpp:48
#define VIENNACL_AMG_MAX_LEVELS
Definition: amg.hpp:54
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
void smooth_jacobi(int level, int const iterations, VectorType &x, VectorType const &rhs) const
(Weighted) Jacobi Smoother (CPU version)
Definition: amg.hpp:526
void clear()
Resets all entries to zero. Does not change the size of the vector.
Definition: vector.hpp:885
void lu_factorize(matrix< SCALARTYPE, viennacl::row_major > &A)
LU factorization of a row-major dense matrix.
Definition: lu.hpp:42
void setup()
Start setup phase for this class and copy data structures.
Definition: amg.hpp:622
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
A tag for algebraic multigrid (AMG). Used to transport information from the user to the implementatio...
Definition: amg_base.hpp:61
void enqueue(KernelType &k, viennacl::ocl::command_queue const &queue)
Enqueues a kernel in the provided queue.
Definition: enqueue.hpp:48
Main kernel class for generating OpenCL kernels for compressed_matrix.
Definition: compressed_matrix.hpp:1039
Implementations of several variants of the AMG coarsening procedure (setup phase). Experimental.
void amg_galerkin_prod(SparseMatrixType &A, SparseMatrixType &P, SparseMatrixType &RES)
Sparse Galerkin product: Calculates RES = trans(P)*A*P.
Definition: amg_base.hpp:1359
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
unsigned int get_coarselevels() const
Definition: amg_base.hpp:111
const_sparse_matrix_adapted_iterator< SCALARTYPE, SizeType,!is_iterator1, true > begin() const
Definition: adapter.hpp:152
static void init(viennacl::ocl::context &ctx)
Definition: compressed_matrix.hpp:1046
A non-const iterator for sparse matrices of type std::vector<std::map<SizeType, SCALARTYPE> > ...
Definition: adapter.hpp:230
AMG preconditioner class, can be supplied to solve()-routines.
Definition: amg.hpp:339
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:692
void lu_substitute(matrix< SCALARTYPE, F1, ALIGNMENT_A > const &A, matrix< SCALARTYPE, F2, ALIGNMENT_B > &B)
LU substitution for the system LU = rhs.
Definition: lu.hpp:201
amg_precond(MatrixType const &mat, amg_tag const &tag)
The constructor. Saves system matrix, tag and builds data structures for setup.
Definition: amg.hpp:376
A class for the sparse matrix type. Uses vector of maps as data structure for higher performance and ...
Definition: amg_base.hpp:376
Implementations of dense direct solvers are found here.
void amg_setup_apply(InternalVectorType &result, InternalVectorType &rhs, InternalVectorType &residual, SparseMatrixType const &A, amg_tag const &tag)
Setup data structures for precondition phase.
Definition: amg.hpp:259
amg_precond()
Definition: amg.hpp:370
void printvector(VectorType const &)
Definition: amg_debug.hpp:80
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
void amg_transform_gpu(InternalType1 &A, InternalType1 &P, InternalType1 &R, InternalType2 &A_setup, InternalType2 &P_setup, amg_tag &tag, viennacl::context ctx)
Save operators after setup phase for GPU computation.
Definition: amg.hpp:219
void init_apply() const
Prepare data structures for preconditioning: Build data structures for precondition phase...
Definition: amg.hpp:636
void apply(VectorType &vec) const
Precondition Operation.
Definition: amg.hpp:445
amg_tag & tag()
Definition: amg.hpp:557
void smooth_jacobi(int level, unsigned int iterations, VectorType &x, VectorType const &rhs) const
Jacobi Smoother (GPU version)
Definition: amg.hpp:768
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:428
Helper classes and functions for the AMG preconditioner. Experimental.
void amg_coarse(unsigned int level, InternalType1 &A, InternalType2 &Pointvector, InternalType3 &Slicing, amg_tag &tag)
Calls the right coarsening procedure.
Definition: amg_coarse.hpp:52
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
ScalarType calc_complexity(VectorType &avgstencil)
Returns complexity measures.
Definition: amg.hpp:652
void setup()
Start setup phase for this class and copy data structures.
Definition: amg.hpp:387
A class for the matrix slicing for parallel coarsening schemes (RS0/RS3).
Definition: amg_base.hpp:1168
unsigned int get_presmooth() const
Definition: amg_base.hpp:105
A class for the AMG points. Holds pointers of type amg_point in a vector that can be accessed using [...
Definition: amg_base.hpp:935
Implementations of several variants of the AMG interpolation operators (setup phase). Experimental.
detail::amg::amg_tag amg_tag
Definition: amg.hpp:60
void switch_memory_context(T &obj, viennacl::context new_ctx)
Generic convenience routine for migrating data of an object to a new memory domain.
Definition: memory.hpp:624