ViennaCL - The Vienna Computing Library  1.5.2
amg_base.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_BASE_HPP_
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_BASE_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/operation.hpp>
28 #include <boost/numeric/ublas/vector.hpp>
29 #include <cmath>
30 #include <set>
31 #include <list>
32 #include <algorithm>
33 
34 #include <map>
35 #ifdef VIENNACL_WITH_OPENMP
36 #include <omp.h>
37 #endif
38 
39 #include "amg_debug.hpp"
40 
41 #define VIENNACL_AMG_COARSE_RS 1
42 #define VIENNACL_AMG_COARSE_ONEPASS 2
43 #define VIENNACL_AMG_COARSE_RS0 3
44 #define VIENNACL_AMG_COARSE_RS3 4
45 #define VIENNACL_AMG_COARSE_AG 5
46 #define VIENNACL_AMG_INTERPOL_DIRECT 1
47 #define VIENNACL_AMG_INTERPOL_CLASSIC 2
48 #define VIENNACL_AMG_INTERPOL_AG 3
49 #define VIENNACL_AMG_INTERPOL_SA 4
50 
51 namespace viennacl
52 {
53  namespace linalg
54  {
55  namespace detail
56  {
57  namespace amg
58  {
61  class amg_tag
62  {
63  public:
76  amg_tag(unsigned int coarse = 1,
77  unsigned int interpol = 1,
78  double threshold = 0.25,
79  double interpolweight = 0.2,
80  double jacobiweight = 1,
81  unsigned int presmooth = 1,
82  unsigned int postsmooth = 1,
83  unsigned int coarselevels = 0)
84  : coarse_(coarse), interpol_(interpol),
85  threshold_(threshold), interpolweight_(interpolweight), jacobiweight_(jacobiweight),
86  presmooth_(presmooth), postsmooth_(postsmooth), coarselevels_(coarselevels) {}
87 
88  // Getter-/Setter-Functions
89  void set_coarse(unsigned int coarse) { if (coarse > 0) coarse_ = coarse; }
90  unsigned int get_coarse() const { return coarse_; }
91 
92  void set_interpol(unsigned int interpol) { if (interpol > 0) interpol_ = interpol; }
93  unsigned int get_interpol() const { return interpol_; }
94 
95  void set_threshold(double threshold) { if (threshold > 0 && threshold <= 1) threshold_ = threshold; }
96  double get_threshold() const{ return threshold_; }
97 
98  void set_as(double jacobiweight) { if (jacobiweight > 0 && jacobiweight <= 2) jacobiweight_ = jacobiweight; }
99  double get_interpolweight() const { return interpolweight_; }
100 
101  void set_interpolweight(double interpolweight) { if (interpolweight > 0 && interpolweight <= 2) interpolweight_ = interpolweight; }
102  double get_jacobiweight() const { return jacobiweight_; }
103 
104  void set_presmooth(int presmooth) { if (presmooth >= 0) presmooth_ = presmooth; }
105  unsigned int get_presmooth() const { return presmooth_; }
106 
107  void set_postsmooth(int postsmooth) { if (postsmooth >= 0) postsmooth_ = postsmooth; }
108  unsigned int get_postsmooth() const { return postsmooth_; }
109 
110  void set_coarselevels(int coarselevels) { if (coarselevels >= 0) coarselevels_ = coarselevels; }
111  unsigned int get_coarselevels() const { return coarselevels_; }
112 
113  private:
114  unsigned int coarse_, interpol_;
115  double threshold_, interpolweight_, jacobiweight_;
116  unsigned int presmooth_, postsmooth_, coarselevels_;
117  };
118 
123  template <typename InternalType, typename IteratorType, typename ScalarType>
125  {
126  private:
127  InternalType *m_;
128  IteratorType iter_;
129  unsigned int i_,j_;
130  ScalarType s_;
131 
132  public:
134 
142  amg_nonzero_scalar(InternalType *m,
143  IteratorType & iter,
144  unsigned int i,
145  unsigned int j,
146  ScalarType s = 0): m_(m), iter_(iter), i_(i), j_(j), s_(s) {}
147 
151  ScalarType operator = (const ScalarType value)
152  {
153  s_ = value;
154  // Only write if scalar is nonzero
155  if (s_ == 0) return s_;
156  // Write to m_ using iterator iter_ or indices (i_,j_)
157  m_->addscalar (iter_,i_,j_,s_);
158  return s_;
159  }
160 
164  ScalarType operator += (const ScalarType value)
165  {
166  // If zero is added, then no change necessary
167  if (value == 0)
168  return s_;
169 
170  s_ += value;
171  // Remove entry if resulting scalar is zero
172  if (s_ == 0)
173  {
174  m_->removescalar(iter_,i_);
175  return s_;
176  }
177  //Write to m_ using iterator iter_ or indices (i_,j_)
178  m_->addscalar (iter_,i_,j_,s_);
179  return s_;
180  }
181  ScalarType operator ++ (int)
182  {
183  s_++;
184  if (s_ == 0)
185  m_->removescalar(iter_,i_);
186  m_->addscalar (iter_,i_,j_,s_);
187  return s_;
188  }
189  ScalarType operator ++ ()
190  {
191  s_++;
192  if (s_ == 0)
193  m_->removescalar(iter_,i_);
194  m_->addscalar (iter_,i_,j_,s_);
195  return s_;
196  }
197  operator ScalarType (void) { return s_; }
198  };
199 
202  template <typename InternalType>
204  {
205  private:
207  typedef typename InternalType::mapped_type ScalarType;
208 
209  InternalType & internal_vec;
210  typename InternalType::iterator iter;
211 
212  public:
213 
218  amg_sparsevector_iterator(InternalType & vec, bool begin=true): internal_vec(vec)
219  {
220  if (begin)
221  iter = internal_vec.begin();
222  else
223  iter = internal_vec.end();
224  }
225 
226  bool operator == (self_type other)
227  {
228  if (iter == other.iter)
229  return true;
230  else
231  return false;
232  }
233  bool operator != (self_type other)
234  {
235  if (iter != other.iter)
236  return true;
237  else
238  return false;
239  }
240 
241  self_type & operator ++ () const { iter++; return *this; }
242  self_type & operator ++ () { iter++; return *this; }
243  self_type & operator -- () const { iter--; return *this; }
244  self_type & operator -- () { iter--; return *this; }
245  ScalarType & operator * () const { return (*iter).second; }
246  ScalarType & operator * () { return (*iter).second; }
247  unsigned int index() const { return (*iter).first; }
248  unsigned int index() { return (*iter).first; }
249  };
250 
253  template <typename ScalarType>
255  {
256  public:
257  typedef ScalarType value_type;
258 
259  private:
260  // A map is used internally which saves all non-zero elements with pairs of (index,value)
261  typedef std::map<unsigned int,ScalarType> InternalType;
264 
265  // Size is only a dummy variable. Not needed for internal map structure but for compatible vector interface.
266  unsigned int size_;
267  InternalType internal_vector;
268 
269  public:
271  typedef typename InternalType::const_iterator const_iterator;
272 
273  public:
277  amg_sparsevector(unsigned int size = 0): size_(size)
278  {
279  internal_vector = InternalType();
280  }
281 
282  void resize(unsigned int size) { size_ = size; }
283  unsigned int size() const { return size_;}
284 
285  // Returns number of non-zero entries in vector equal to the size of the underlying map.
286  unsigned int internal_size() const { return static_cast<unsigned int>(internal_vector.size()); }
287  // Delete underlying map.
288  void clear() { internal_vector.clear(); }
289  // Remove entry at position i.
290  void remove(unsigned int i) { internal_vector.erase(i); }
291 
292  // Add s to the entry at position i
293  void add (unsigned int i, ScalarType s)
294  {
295  typename InternalType::iterator iter = internal_vector.find(i);
296  // If there is no entry at position i, add new entry at that position
297  if (iter == internal_vector.end())
298  addscalar(iter,i,i,s);
299  else
300  {
301  s += (*iter).second;
302  // If new value is zero, then erase the entry, otherwise write new value
303  if (s != 0)
304  (*iter).second = s;
305  else
306  internal_vector.erase(iter);
307  }
308  }
309 
310  // Write to the map. Is called from non-zero scalar type.
311  template <typename IteratorType>
312  void addscalar(IteratorType & iter, unsigned int i, unsigned int /* j */, ScalarType s)
313  {
314  // Don't write if value is zero
315  if (s == 0)
316  return;
317 
318  // If entry is already present, overwrite value, otherwise make new entry
319  if (iter != internal_vector.end())
320  (*iter).second = s;
321  else
322  internal_vector[i] = s;
323  }
324 
325  // Remove value from the map. Is called from non-zero scalar type.
326  template <typename IteratorType>
327  void removescalar(IteratorType & iter, unsigned int /* i */) { internal_vector.erase(iter); }
328 
329  // Bracket operator. Returns non-zero scalar type with actual values of the respective entry which calls addscalar/removescalar after value is altered.
330  NonzeroScalarType operator [] (unsigned int i)
331  {
332  typename InternalType::iterator it = internal_vector.find(i);
333  // If value is already present then build non-zero scalar with actual value, otherwise 0.
334  if (it != internal_vector.end())
335  return NonzeroScalarType (this,it,i,i,(*it).second);
336  else
337  return NonzeroScalarType (this,it,i,i,0);
338  }
339 
340  // Use internal data structure directly for read-only access. No need to use non-zero scalar as no write access possible.
341  ScalarType operator [] (unsigned int i) const
342  {
343  const_iterator it = internal_vector.find(i);
344 
345  if (it != internal_vector.end())
346  return (*it).second;
347  else
348  return 0;
349  }
350 
351  // Iterator functions.
352  iterator begin() { return iterator(internal_vector); }
353  const_iterator begin() const { return internal_vector.begin(); }
354  iterator end() { return iterator(internal_vector,false); }
355  const_iterator end() const { return internal_vector.end(); }
356 
357  // checks whether value at index i is nonzero. More efficient than doing [] == 0.
358  bool isnonzero(unsigned int i) const { return internal_vector.find(i) != internal_vector.end(); }
359 
360  // Copies data into a ublas vector type.
361  operator boost::numeric::ublas::vector<ScalarType> (void)
362  {
363  boost::numeric::ublas::vector<ScalarType> vec (size_);
364  for (iterator iter = begin(); iter != end(); ++iter)
365  vec [iter.index()] = *iter;
366  return vec;
367  }
368  };
369 
375  template <typename ScalarType>
377  {
378  public:
379  typedef ScalarType value_type;
380  private:
381  typedef std::map<unsigned int,ScalarType> RowType;
382  typedef std::vector<RowType> InternalType;
384 
385  // Adapter is used for certain functionality, especially iterators.
388 
389  // Non-zero scalar is used to write to the matrix.
391 
392  // Holds matrix coefficients.
393  InternalType internal_mat;
394  // Holds matrix coefficient of transposed matrix if built.
395  // Note: Only internal_mat is written using operators and methods while internal_mat_trans is built from internal_mat using do_trans().
396  InternalType internal_mat_trans;
397  // Saves sizes.
398  vcl_size_t s1, s2;
399 
400  // True if the transposed of the matrix is used (for calculations, iteration, etc.).
401  bool transposed_mode;
402  // True if the transposed is already built (saved in internal_mat_trans) and also up to date (no changes to internal_mat).
403  bool transposed;
404 
405  public:
410 
413  {
414  transposed_mode = false;
415  transposed = false;
416  }
417 
422  amg_sparsematrix (unsigned int i, unsigned int j)
423  {
424  AdapterType a (internal_mat, i, j);
425  a.resize (i,j,false);
426  AdapterType a_trans (internal_mat_trans, j, i);
427  a_trans.resize (j,i,false);
428  s1 = i;
429  s2 = j;
430  a.clear();
431  a_trans.clear();
432  transposed_mode = false;
433  transposed = false;
434  }
435 
440  amg_sparsematrix (std::vector<std::map<unsigned int, ScalarType> > const & mat)
441  {
442  AdapterType a (internal_mat, mat.size(), mat.size());
443  AdapterType a_trans (internal_mat_trans, mat.size(), mat.size());
444  a.resize(mat.size(), mat.size());
445  a_trans.resize(mat.size(), mat.size());
446 
447  internal_mat = mat;
448  s1 = s2 = mat.size();
449 
450  transposed_mode = false;
451  transposed = false;
452  }
453 
458  template <typename MatrixType>
459  amg_sparsematrix (MatrixType const & mat)
460  {
461  AdapterType a (internal_mat, mat.size1(), mat.size2());
462  AdapterType a_trans (internal_mat_trans, mat.size2(), mat.size1());
463  a.resize(mat.size1(), mat.size2());
464  a_trans.resize(mat.size2(), mat.size1());
465  s1 = mat.size1();
466  s2 = mat.size2();
467  a.clear();
468  a_trans.clear();
469 
470  for (typename MatrixType::const_iterator1 row_iter = mat.begin1(); row_iter != mat.end1(); ++row_iter)
471  {
472  for (typename MatrixType::const_iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
473  {
474  if (*col_iter != 0)
475  {
476  unsigned int x = static_cast<unsigned int>(col_iter.index1());
477  unsigned int y = static_cast<unsigned int>(col_iter.index2());
478  a (x,y) = *col_iter;
479  a_trans (y,x) = *col_iter;
480  }
481  }
482  }
483  transposed_mode = false;
484  transposed = true;
485  }
486 
487  // Build transposed of the current matrix.
488  void do_trans()
489  {
490  // Do it only once if called in a parallel section
491  #ifdef VIENNACL_WITH_OPENMP
492  #pragma omp critical
493  #endif
494  {
495  // Only build transposed if it is not built or not up to date
496  if (!transposed)
497  {
498  // Mode has to be set to standard mode temporarily
499  bool save_mode = transposed_mode;
500  transposed_mode = false;
501 
502  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
503  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
504  internal_mat_trans[col_iter.index2()][static_cast<unsigned int>(col_iter.index1())] = *col_iter;
505 
506  transposed_mode = save_mode;
507  transposed = true;
508  }
509  }
510  } //do_trans()
511 
512  // Set transposed mode (true=transposed, false=regular)
513  void set_trans(bool mode)
514  {
515  transposed_mode = mode;
516  if (mode)
517  do_trans();
518  }
519 
520  bool get_trans() const { return transposed_mode; }
521 
522  // Checks whether coefficient (i,j) is non-zero. More efficient than using (i,j) == 0.
523  bool isnonzero (unsigned int i, unsigned int j) const
524  {
525  if (!transposed_mode)
526  {
527  if (internal_mat[i].find(j) != internal_mat[i].end())
528  return true;
529  else
530  return false;
531  }
532  else
533  {
534  if (internal_mat[j].find(i) != internal_mat[j].end())
535  return true;
536  else
537  return false;
538  }
539  } //isnonzero()
540 
541  // Add s to value at (i,j)
542  void add (unsigned int i, unsigned int j, ScalarType s)
543  {
544  // If zero is added then do nothing.
545  if (s == 0)
546  return;
547 
548  typename RowType::iterator col_iter = internal_mat[i].find(j);
549  // If there is no entry at position (i,j), then make new entry.
550  if (col_iter == internal_mat[i].end())
551  addscalar(col_iter,i,j,s);
552  else
553  {
554  s += (*col_iter).second;
555  // Update value and erase entry if value is zero.
556  if (s != 0)
557  (*col_iter).second = s;
558  else
559  internal_mat[i].erase(col_iter);
560  }
561  transposed = false;
562  } //add()
563 
564  // Write to the internal data structure. Is called from non-zero scalar type.
565  template <typename IteratorType>
566  void addscalar(IteratorType & iter, unsigned int i, unsigned int j, ScalarType s)
567  {
568  // Don't write if value is zero
569  if (s == 0)
570  return;
571 
572  if (iter != internal_mat[i].end())
573  (*iter).second = s;
574  else
575  internal_mat[i][j] = s;
576 
577  transposed = false;
578  }
579 
580  // Remove entry from internal data structure. Is called from non-zero scalar type.
581  template <typename IteratorType>
582  void removescalar(IteratorType & iter, unsigned int i)
583  {
584  internal_mat[i].erase(iter);
585  transposed = false;
586  }
587 
588  // Return non-zero scalar at position (i,j). Value is written to the non-zero scalar and updated via addscalar()/removescalar().
589  NonzeroScalarType operator()(unsigned int i, unsigned int j)
590  {
591  typename RowType::iterator iter;
592 
593  if (!transposed_mode)
594  {
595  iter = internal_mat[i].find(j);
596  if (iter != internal_mat[i].end())
597  return NonzeroScalarType (this,iter,i,j,(*iter).second);
598  else
599  return NonzeroScalarType (this,iter,i,j,0);
600  }
601  else
602  {
603  iter = internal_mat[j].find(i);
604  if (iter != internal_mat[j].end())
605  return NonzeroScalarType (this,iter,j,i,(*iter).second);
606  else
607  return NonzeroScalarType (this,iter,j,i,0);
608  }
609  }
610 
611  // For read-only access return the actual value directly. Non-zero datatype not needed as no write access possible.
612  ScalarType operator()(unsigned int i, unsigned int j) const
613  {
614  typename RowType::const_iterator iter;
615 
616  if (!transposed_mode)
617  {
618  iter = internal_mat[i].find(j);
619  if (iter != internal_mat[i].end())
620  return (*iter).second;
621  else
622  return 0;
623  }
624  else
625  {
626  iter = internal_mat[j].find(i);
627  if (iter != internal_mat[j].end())
628  return (*iter).second;
629  else
630  return 0;
631  }
632  }
633 
634  void resize(unsigned int i, unsigned int j, bool preserve = true)
635  {
636  AdapterType a (internal_mat);
637  a.resize(i,j,preserve);
638  AdapterType a_trans (internal_mat_trans);
639  a_trans.resize(j,i,preserve);
640  s1 = i;
641  s2 = j;
642  }
643 
644  void clear()
645  {
646  AdapterType a (internal_mat, s1, s2);
647  a.clear();
648  AdapterType a_trans (internal_mat_trans, s2, s1);
649  a_trans.clear();
650  transposed = true;
651  }
652 
654  {
655  if (!transposed_mode)
656  return s1;
657  else
658  return s2;
659  }
660 
662  {
663  if (!transposed_mode)
664  return s1;
665  else
666  return s2;
667  }
668 
669 
671  {
672  if (!transposed_mode)
673  return s2;
674  else
675  return s1;
676  }
677 
679  {
680  if (!transposed_mode)
681  return s2;
682  else
683  return s1;
684  }
685 
686  iterator1 begin1(bool trans = false)
687  {
688  if (!trans && !transposed_mode)
689  {
690  AdapterType a (internal_mat, s1, s2);
691  return a.begin1();
692  }
693  else
694  {
695  do_trans();
696  AdapterType a_trans (internal_mat_trans, s2, s1);
697  return a_trans.begin1();
698  }
699  }
700 
701  iterator1 end1(bool trans = false)
702  {
703  if (!trans && !transposed_mode)
704  {
705  AdapterType a (internal_mat, s1, s2);
706  return a.end1();
707  }
708  else
709  {
710  //do_trans();
711  AdapterType a_trans (internal_mat_trans, s2, s1);
712  return a_trans.end1();
713  }
714  }
715 
716  iterator2 begin2(bool trans = false)
717  {
718  if (!trans && !transposed_mode)
719  {
720  AdapterType a (internal_mat, s1, s2);
721  return a.begin2();
722  }
723  else
724  {
725  do_trans();
726  AdapterType a_trans (internal_mat_trans, s2, s1);
727  return a_trans.begin2();
728  }
729  }
730 
731  iterator2 end2(bool trans = false)
732  {
733  if (!trans && !transposed_mode)
734  {
735  AdapterType a (internal_mat, s1, s2);
736  return a.end2();
737  }
738  else
739  {
740  //do_trans();
741  AdapterType a_trans (internal_mat_trans, s2, s1);
742  return a_trans.end2();
743  }
744  }
745 
746  const_iterator1 begin1() const
747  {
748  // Const_iterator of transposed can only be used if transposed matrix is already built and up to date.
749  assert((!transposed_mode || (transposed_mode && transposed)) && bool("Error: Cannot build const_iterator when transposed has not been built yet!"));
750  ConstAdapterType a_const (internal_mat, s1, s2);
751  return a_const.begin1();
752  }
753 
754  const_iterator1 end1(bool trans = false) const
755  {
756  assert((!transposed_mode || (transposed_mode && transposed)) && bool("Error: Cannot build const_iterator when transposed has not been built yet!"));
757  ConstAdapterType a_const (internal_mat, trans ? s2 : s1, trans ? s1 : s2);
758  return a_const.end1();
759  }
760 
761  const_iterator2 begin2(bool trans = false) const
762  {
763  assert((!transposed_mode || (transposed_mode && transposed)) && bool("Error: Cannot build const_iterator when transposed has not been built yet!"));
764  ConstAdapterType a_const (internal_mat, trans ? s2 : s1, trans ? s1 : s2);
765  return a_const.begin2();
766  }
767 
768  const_iterator2 end2(bool trans = false) const
769  {
770  assert((!transposed_mode || (transposed_mode && transposed)) && bool("Error: Cannot build const_iterator when transposed has not been built yet!"));
771  ConstAdapterType a_const (internal_mat, trans ? s2 : s1, trans ? s1 : s2);
772  return a_const.end2();
773  }
774 
775  // Returns pointer to the internal data structure. Improves performance of copy operation to GPU.
776  std::vector<std::map<unsigned int, ScalarType> > * get_internal_pointer()
777  {
778  if (!transposed_mode)
779  return &internal_mat;
780 
781  if (!transposed)
782  do_trans();
783  return &internal_mat_trans;
784  }
785  operator boost::numeric::ublas::compressed_matrix<ScalarType> (void)
786  {
787  boost::numeric::ublas::compressed_matrix<ScalarType> mat;
788  mat.resize(size1(),size2(),false);
789  mat.clear();
790 
791  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
792  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
793  mat (col_iter.index1(), col_iter.index2()) = *col_iter;
794 
795  return mat;
796  }
797  operator boost::numeric::ublas::matrix<ScalarType> (void)
798  {
799  boost::numeric::ublas::matrix<ScalarType> mat;
800  mat.resize(size1(),size2(),false);
801  mat.clear();
802 
803  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
804  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
805  mat (col_iter.index1(), col_iter.index2()) = *col_iter;
806 
807  return mat;
808  }
809  };
810 
816  class amg_point
817  {
818  private:
820 
821  unsigned int index_;
822  unsigned int influence_;
823  // Determines whether point is undecided.
824  bool undecided_;
825  // Determines wheter point is C point (true) or F point (false).
826  bool cpoint_;
827  unsigned int coarse_index_;
828  // Index offset of parallel coarsening. In that case a point acts as if it had an index of index_-offset_ and treats other points as if they had an index of index+offset_
829  unsigned int offset_;
830  // Aggregate the point belongs to.
831  unsigned int aggregate_;
832 
833  // Holds all points influencing this point.
834  ListType influencing_points;
835  // Holds all points that are influenced by this point.
836  ListType influenced_points;
837 
838  public:
841 
844  amg_point (unsigned int index, unsigned int size): index_(index), influence_(0), undecided_(true), cpoint_(false), coarse_index_(0), offset_(0), aggregate_(0)
845  {
846  influencing_points = ListType(size);
847  influenced_points = ListType(size);
848  }
849 
850  void set_offset(unsigned int offset) { offset_ = offset; }
851  unsigned int get_offset() { return offset_; }
852  void set_index(unsigned int index) { index_ = index+offset_; }
853  unsigned int get_index() const { return index_-offset_; }
854  unsigned int get_influence() const { return influence_; }
855  void set_aggregate(unsigned int aggregate) { aggregate_ = aggregate; }
856  unsigned int get_aggregate () { return aggregate_; }
857 
858  bool is_cpoint() const { return cpoint_ && !undecided_; }
859  bool is_fpoint() const { return !cpoint_ && !undecided_; }
860  bool is_undecided() const { return undecided_; }
861 
862  // Returns number of influencing points
863  unsigned int number_influencing() const { return influencing_points.internal_size(); }
864  // Returns true if *point is influencing this point
865  bool is_influencing(amg_point* point) const { return influencing_points.isnonzero(point->get_index()+offset_); }
866  // Add *point to influencing points
867  void add_influencing_point(amg_point* point) { influencing_points[point->get_index()+offset_] = point; }
868  // Add *point to influenced points
869  void add_influenced_point(amg_point* point) { influenced_points[point->get_index()+offset_] = point; }
870 
871  // Clear influencing points
872  void clear_influencing() { influencing_points.clear(); }
873  // Clear influenced points
874  void clear_influenced() {influenced_points.clear(); }
875 
876 
877  unsigned int get_coarse_index() const { return coarse_index_; }
878  void set_coarse_index(unsigned int index) { coarse_index_ = index; }
879 
880  // Calculates the initial influence measure equal to the number of influenced points.
881  void calc_influence() { influence_ = influenced_points.internal_size(); }
882 
883  // Add to influence measure.
884  unsigned int add_influence(int add)
885  {
886  influence_ += add;
887  return influence_;
888  }
889  // Make this point C point. Only call via amg_pointvector.
890  void make_cpoint()
891  {
892  undecided_ = false;
893  cpoint_ = true;
894  influence_ = 0;
895  }
896  // Make this point F point. Only call via amg_pointvector.
897  void make_fpoint()
898  {
899  undecided_ = false;
900  cpoint_ = false;
901  influence_ = 0;
902  }
903  // Switch point from F to C point. Only call via amg_pointvector.
904  void switch_ftoc() { cpoint_ = true; }
905 
906  // Iterator handling for influencing and influenced points.
907  iterator begin_influencing() { return influencing_points.begin(); }
908  iterator end_influencing() { return influencing_points.end(); }
909  const_iterator begin_influencing() const { return influencing_points.begin(); }
910  const_iterator end_influencing() const { return influencing_points.end(); }
911  iterator begin_influenced() { return influenced_points.begin(); }
912  iterator end_influenced() { return influenced_points.end(); }
913  const_iterator begin_influenced() const { return influenced_points.begin(); }
914  const_iterator end_influenced() const { return influenced_points.end(); }
915  };
916 
919  struct classcomp
920  {
921  // Function returns true if l comes before r in the ordering.
922  bool operator() (amg_point* l, amg_point* r) const
923  {
924  // Map is sorted by influence number starting with the highest
925  // If influence number is the same then lowest point index comes first
926  return (l->get_influence() < r->get_influence() || (l->get_influence() == r->get_influence() && l->get_index() > r->get_index()));
927  }
928  };
929 
936  {
937  private:
938  // Type for the sorted list
939  typedef std::set<amg_point*,classcomp> ListType;
940  // Type for the vector of pointers
941  typedef std::vector<amg_point*> VectorType;
942  VectorType pointvector;
943  ListType pointlist;
944  unsigned int size_;
945  unsigned int c_points, f_points;
946 
947  public:
948  typedef VectorType::iterator iterator;
949  typedef VectorType::const_iterator const_iterator;
950 
954  amg_pointvector(unsigned int size = 0): size_(size)
955  {
956  pointvector = VectorType(size);
957  c_points = f_points = 0;
958  }
959 
960  // Construct all the points dynamically and save pointers into vector.
961  void init_points()
962  {
963  for (unsigned int i=0; i<size(); ++i)
964  pointvector[i] = new amg_point(i,size());
965  }
966  // Delete all the points.
968  {
969  for (unsigned int i=0; i<size(); ++i)
970  delete pointvector[i];
971  }
972  // Add point to the vector. Note: User has to make sure that point at point->get_index() does not exist yet, otherwise it will be overwritten!
973  void add_point(amg_point *point)
974  {
975  pointvector[point->get_index()] = point;
976  if (point->is_cpoint()) c_points++;
977  else if (point->is_fpoint()) f_points++;
978  }
979 
980  // Update C and F count for point *point.
981  // Necessary if C and F points were constructed outside this data structure (e.g. by parallel coarsening RS0 or RS3).
982  void update_cf(amg_point *point)
983  {
984  if (point->is_cpoint()) c_points++;
985  else if (point->is_fpoint()) f_points++;
986  }
987  // Clear the C and F point count.
988  void clear_cf() { c_points = f_points = 0; }
989 
990  // Clear both point lists.
992  {
993  for (iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
994  {
995  (*iter)->clear_influencing();
996  (*iter)->clear_influenced();
997  }
998  }
999 
1000  amg_point* operator [] (unsigned int i) const { return pointvector[i]; }
1001  iterator begin() { return pointvector.begin(); }
1002  iterator end() { return pointvector.end(); }
1003  const_iterator begin() const { return pointvector.begin(); }
1004  const_iterator end() const { return pointvector.end(); }
1005 
1006  void resize(unsigned int size)
1007  {
1008  size_ = size;
1009  pointvector = VectorType(size);
1010  }
1011  unsigned int size() const { return size_; }
1012 
1013  // Returns number of C points
1014  unsigned int get_cpoints() const { return c_points; }
1015  // Returns number of F points
1016  unsigned int get_fpoints() const { return f_points; }
1017 
1018  // Does the initial sorting of points into the list. Sorting is automatically done by the std::set data type.
1019  void sort()
1020  {
1021  for (iterator iter = begin(); iter != end(); ++iter)
1022  pointlist.insert(*iter);
1023  }
1024  // Returns the point with the highest influence measure
1026  {
1027  // No points remaining? Return NULL such that coarsening will stop.
1028  if (pointlist.size() == 0)
1029  return NULL;
1030  // If point with highest influence measure (end of the list) has measure of zero, then no further C points can be constructed. Return NULL.
1031  if ((*(--pointlist.end()))->get_influence() == 0)
1032  return NULL;
1033  // Otherwise, return the point with highest influence measure located at the end of the list.
1034  else
1035  return *(--pointlist.end());
1036  }
1037  // Add "add" to influence measure for point *point in the sorted list.
1038  void add_influence(amg_point* point, unsigned int add)
1039  {
1040  ListType::iterator iter = pointlist.find(point);
1041  // If point is not in the list then stop.
1042  if (iter == pointlist.end()) return;
1043 
1044  // Save iterator and decrement
1045  ListType::iterator iter2 = iter;
1046  iter2--;
1047 
1048  // Point has to be erased first as changing the value does not re-order the std::set
1049  pointlist.erase(iter);
1050  point->add_influence(add);
1051 
1052  // Insert point back into the list. Using the iterator improves performance. The new position has to be at the same position or to the right of the old.
1053  pointlist.insert(iter2,point);
1054  }
1055  // Make *point to C point and remove from sorted list
1056  void make_cpoint(amg_point* point)
1057  {
1058  pointlist.erase(point);
1059  point->make_cpoint();
1060  c_points++;
1061  }
1062  // Make *point to F point and remove from sorted list
1063  void make_fpoint(amg_point* point)
1064  {
1065  pointlist.erase(point);
1066  point->make_fpoint();
1067  f_points++;
1068  }
1069  // Swich *point from F to C point
1070  void switch_ftoc(amg_point* point)
1071  {
1072  point->switch_ftoc();
1073  c_points++;
1074  f_points--;
1075  }
1076 
1077  // Build vector of indices for C point on the coarse level.
1079  {
1080  unsigned int count = 0;
1081  // Use simple counter for index creation.
1082  for (iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
1083  {
1084  // Set index on coarse level using counter variable
1085  if ((*iter)->is_cpoint())
1086  {
1087  (*iter)->set_coarse_index(count);
1088  count++;
1089  }
1090  }
1091  }
1092 
1093  // Return information for debugging purposes
1094  template <typename MatrixType>
1095  void get_influence_matrix(MatrixType & mat) const
1096  {
1097  mat = MatrixType(size(),size());
1098  mat.clear();
1099 
1100  for (const_iterator row_iter = begin(); row_iter != end(); ++row_iter)
1101  for (amg_point::iterator col_iter = (*row_iter)->begin_influencing(); col_iter != (*row_iter)->end_influencing(); ++col_iter)
1102  mat((*row_iter)->get_index(),(*col_iter)->get_index()) = true;
1103  }
1104  template <typename VectorType>
1105  void get_influence(VectorType & vec) const
1106  {
1107  vec = VectorType(size_);
1108  vec.clear();
1109 
1110  for (const_iterator iter = begin(); iter != end(); ++iter)
1111  vec[(*iter)->get_index()] = (*iter)->get_influence();
1112  }
1113  template <typename VectorType>
1114  void get_sorting(VectorType & vec) const
1115  {
1116  vec = VectorType(pointlist.size());
1117  vec.clear();
1118  unsigned int i=0;
1119 
1120  for (ListType::const_iterator iter = pointlist.begin(); iter != pointlist.end(); ++iter)
1121  {
1122  vec[i] = (*iter)->get_index();
1123  i++;
1124  }
1125  }
1126  template <typename VectorType>
1127  void get_C(VectorType & vec) const
1128  {
1129  vec = VectorType(size_);
1130  vec.clear();
1131 
1132  for (const_iterator iter = begin(); iter != end(); ++iter)
1133  {
1134  if ((*iter)->is_cpoint())
1135  vec[(*iter)->get_index()] = true;
1136  }
1137  }
1138  template <typename VectorType>
1139  void get_F(VectorType & vec) const
1140  {
1141  vec = VectorType(size_);
1142  vec.clear();
1143 
1144  for (const_iterator iter = begin(); iter != end(); ++iter)
1145  {
1146  if ((*iter)->is_fpoint())
1147  vec[(*iter)->get_index()] = true;
1148  }
1149  }
1150  template <typename MatrixType>
1151  void get_Aggregates(MatrixType & mat) const
1152  {
1153  mat = MatrixType(size_,size_);
1154  mat.clear();
1155 
1156  for (const_iterator iter = begin(); iter != end(); ++iter)
1157  {
1158  if (!(*iter)->is_undecided())
1159  mat((*iter)->get_aggregate(),(*iter)->get_index()) = true;
1160  }
1161  }
1162  };
1163 
1167  template <typename InternalType1, typename InternalType2>
1169  {
1170  typedef typename InternalType1::value_type SparseMatrixType;
1171  typedef typename InternalType2::value_type PointVectorType;
1172 
1173  public:
1174  // Data structures on a per-processor basis.
1175  boost::numeric::ublas::vector<InternalType1> A_slice;
1176  boost::numeric::ublas::vector<InternalType2> Pointvector_slice;
1177  // Holds the offsets showing the indices for which a new slice begins.
1178  boost::numeric::ublas::vector<boost::numeric::ublas::vector<unsigned int> > Offset;
1179 
1180  unsigned int threads_;
1181  unsigned int levels_;
1182 
1183  void init(unsigned int levels, unsigned int threads = 0)
1184  {
1185  // Either use the number of threads chosen by the user or the maximum number of threads available on the processor.
1186  if (threads == 0)
1187  #ifdef VIENNACL_WITH_OPENMP
1188  threads_ = omp_get_num_procs();
1189  #else
1190  threads_ = 1;
1191  #endif
1192  else
1193  threads_ = threads;
1194 
1195  levels_ = levels;
1196 
1197  A_slice.resize(threads_);
1198  Pointvector_slice.resize(threads_);
1199  // Offset has threads_+1 entries to also hold the total size
1200  Offset.resize(threads_+1);
1201 
1202  for (unsigned int i=0; i<threads_; ++i)
1203  {
1204  A_slice[i].resize(levels_);
1205  Pointvector_slice[i].resize(levels_);
1206  // Offset needs one more level for the build-up of the next offset
1207  Offset[i].resize(levels_+1);
1208  }
1209  Offset[threads_].resize(levels_+1);
1210  } //init()
1211 
1212  // Slice matrix A into as many parts as threads are used.
1213  void slice (unsigned int level, InternalType1 const & A, InternalType2 const & Pointvector)
1214  {
1215  // On the finest level, build a new slicing first.
1216  if (level == 0)
1217  slice_new (level, A);
1218 
1219  // On coarser levels use the same slicing as on the finest level (Points stay together on the same thread on all levels).
1220  // This is necessary as due to interpolation and galerkin product there only exist connections between points on the same thread on coarser levels.
1221  // Note: Offset is determined in amg_coarse_rs0() after fine level was built.
1222  slice_build (level, A, Pointvector);
1223  }
1224 
1225  // Join point data structure into Pointvector
1226  void join (unsigned int level, InternalType2 & Pointvector) const
1227  {
1228  typedef typename InternalType2::value_type PointVectorType;
1229 
1230  // Reset index offset of all points and update overall C and F point count
1231  Pointvector[level].clear_cf();
1232  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
1233  {
1234  (*iter)->set_offset(0);
1235  Pointvector[level].update_cf(*iter);
1236  }
1237  }
1238 
1239  private:
1244  void slice_new (unsigned int level, InternalType1 const & A)
1245  {
1246  // Determine index offset of all the slices (index of A[level] when the respective slice starts).
1247  #ifdef VIENNACL_WITH_OPENMP
1248  #pragma omp parallel for
1249  #endif
1250  for (long i=0; i<=static_cast<long>(threads_); ++i)
1251  {
1252  // Offset of first piece is zero. Pieces 1,...,threads-1 have equal size while the last one might be greater.
1253  if (i == 0) Offset[i][level] = 0;
1254  else if (i == threads_) Offset[i][level] = static_cast<unsigned int>(A[level].size1());
1255  else Offset[i][level] = static_cast<unsigned int>(i*(A[level].size1()/threads_));
1256  }
1257  }
1258 
1264  void slice_build (unsigned int level, InternalType1 const & A, InternalType2 const & Pointvector)
1265  {
1266  typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
1267  typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
1268 
1269  unsigned int x, y;
1270  amg_point *point;
1271 
1272  #ifdef VIENNACL_WITH_OPENMP
1273  #pragma omp parallel for private (x,y,point)
1274  #endif
1275  for (long i=0; i<static_cast<long>(threads_); ++i)
1276  {
1277  // Allocate space for the matrix slice and the pointvector.
1278  A_slice[i][level] = SparseMatrixType(Offset[i+1][level]-Offset[i][level],Offset[i+1][level]-Offset[i][level]);
1279  Pointvector_slice[i][level] = PointVectorType(Offset[i+1][level]-Offset[i][level]);
1280 
1281  // Iterate over the part that belongs to thread i (from Offset[i][level] to Offset[i+1][level]).
1282  ConstRowIterator row_iter = A[level].begin1();
1283  row_iter += Offset[i][level];
1284  x = static_cast<unsigned int>(row_iter.index1());
1285 
1286  while (x < Offset[i+1][level] && row_iter != A[level].end1())
1287  {
1288  // Set offset for point index and save point for the respective thread
1289  point = Pointvector[level][x];
1290  point->set_offset(Offset[i][level]);
1291  Pointvector_slice[i][level].add_point(point);
1292 
1293  ConstColIterator col_iter = row_iter.begin();
1294  y = static_cast<unsigned int>(col_iter.index2());
1295 
1296  // Save all coefficients from the matrix slice
1297  while (y < Offset[i+1][level] && col_iter != row_iter.end())
1298  {
1299  if (y >= Offset[i][level])
1300  A_slice[i][level](x-Offset[i][level],y-Offset[i][level]) = *col_iter;
1301 
1302  ++col_iter;
1303  y = static_cast<unsigned int>(col_iter.index2());
1304  }
1305 
1306  ++row_iter;
1307  x = static_cast<unsigned int>(row_iter.index1());
1308  }
1309  }
1310  }
1311  };
1312 
1318  template <typename SparseMatrixType>
1319  void amg_mat_prod (SparseMatrixType & A, SparseMatrixType & B, SparseMatrixType & RES)
1320  {
1321  typedef typename SparseMatrixType::value_type ScalarType;
1322  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
1323  typedef typename SparseMatrixType::iterator2 InternalColIterator;
1324 
1325  long x,y,z;
1326  ScalarType prod;
1327  RES = SparseMatrixType(static_cast<unsigned int>(A.size1()), static_cast<unsigned int>(B.size2()));
1328  RES.clear();
1329 
1330  #ifdef VIENNACL_WITH_OPENMP
1331  #pragma omp parallel for private (x,y,z,prod)
1332  #endif
1333  for (x=0; x<static_cast<long>(A.size1()); ++x)
1334  {
1335  InternalRowIterator row_iter = A.begin1();
1336  row_iter += x;
1337  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
1338  {
1339  y = static_cast<unsigned int>(col_iter.index2());
1340  InternalRowIterator row_iter2 = B.begin1();
1341  row_iter2 += y;
1342 
1343  for(InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
1344  {
1345  z = static_cast<unsigned int>(col_iter2.index2());
1346  prod = *col_iter * *col_iter2;
1347  RES.add(x,z,prod);
1348  }
1349  }
1350  }
1351  }
1352 
1358  template <typename SparseMatrixType>
1359  void amg_galerkin_prod (SparseMatrixType & A, SparseMatrixType & P, SparseMatrixType & RES)
1360  {
1361  typedef typename SparseMatrixType::value_type ScalarType;
1362  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
1363  typedef typename SparseMatrixType::iterator2 InternalColIterator;
1364 
1365  long x,y1,y2,z;
1367  RES = SparseMatrixType(static_cast<unsigned int>(P.size2()), static_cast<unsigned int>(P.size2()));
1368  RES.clear();
1369 
1370  #ifdef VIENNACL_WITH_OPENMP
1371  #pragma omp parallel for private (x,y1,y2,z,row)
1372  #endif
1373  for (x=0; x<static_cast<long>(P.size2()); ++x)
1374  {
1375  row = amg_sparsevector<ScalarType>(static_cast<unsigned int>(A.size2()));
1376  InternalRowIterator row_iter = P.begin1(true);
1377  row_iter += x;
1378 
1379  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
1380  {
1381  y1 = static_cast<long>(col_iter.index2());
1382  InternalRowIterator row_iter2 = A.begin1();
1383  row_iter2 += y1;
1384 
1385  for(InternalColIterator col_iter2 = row_iter2.begin(); col_iter2 != row_iter2.end(); ++col_iter2)
1386  {
1387  y2 = static_cast<long>(col_iter2.index2());
1388  row.add (y2, *col_iter * *col_iter2);
1389  }
1390  }
1391  for (typename amg_sparsevector<ScalarType>::iterator iter = row.begin(); iter != row.end(); ++iter)
1392  {
1393  y2 = iter.index();
1394  InternalRowIterator row_iter3 = P.begin1();
1395  row_iter3 += y2;
1396 
1397  for (InternalColIterator col_iter3 = row_iter3.begin(); col_iter3 != row_iter3.end(); ++col_iter3)
1398  {
1399  z = static_cast<long>(col_iter3.index2());
1400  RES.add (x, z, *col_iter3 * *iter);
1401  }
1402  }
1403  }
1404 
1405  #ifdef VIENNACL_AMG_DEBUG
1406  std::cout << "Galerkin Operator: " << std::endl;
1407  printmatrix (RES);
1408  #endif
1409  }
1410 
1416  template <typename SparseMatrixType>
1417  void test_triplematprod(SparseMatrixType & A, SparseMatrixType & P, SparseMatrixType & A_i1)
1418  {
1419  typedef typename SparseMatrixType::value_type ScalarType;
1420 
1421  boost::numeric::ublas::compressed_matrix<ScalarType> A_temp (A.size1(), A.size2());
1422  A_temp = A;
1423  boost::numeric::ublas::compressed_matrix<ScalarType> P_temp (P.size1(), P.size2());
1424  P_temp = P;
1425  P.set_trans(true);
1426  boost::numeric::ublas::compressed_matrix<ScalarType> R_temp (P.size1(), P.size2());
1427  R_temp = P;
1428  P.set_trans(false);
1429 
1430  boost::numeric::ublas::compressed_matrix<ScalarType> RA (R_temp.size1(),A_temp.size2());
1431  RA = boost::numeric::ublas::prod(R_temp,A_temp);
1432  boost::numeric::ublas::compressed_matrix<ScalarType> RAP (RA.size1(),P_temp.size2());
1433  RAP = boost::numeric::ublas::prod(RA,P_temp);
1434 
1435  for (unsigned int x=0; x<RAP.size1(); ++x)
1436  {
1437  for (unsigned int y=0; y<RAP.size2(); ++y)
1438  {
1439  if (std::fabs(static_cast<ScalarType>(RAP(x,y)) - static_cast<ScalarType>(A_i1(x,y))) > 0.0001)
1440  std::cout << x << " " << y << " " << RAP(x,y) << " " << A_i1(x,y) << std::endl;
1441  }
1442  }
1443  }
1444 
1450  template <typename SparseMatrixType, typename PointVectorType>
1451  void test_interpolation(SparseMatrixType & A, SparseMatrixType & P, PointVectorType & Pointvector)
1452  {
1453  for (unsigned int i=0; i<P.size1(); ++i)
1454  {
1455  if (Pointvector.is_cpoint(i))
1456  {
1457  bool set = false;
1458  for (unsigned int j=0; j<P.size2(); ++j)
1459  {
1460  if (P.isnonzero(i,j))
1461  {
1462  if (P(i,j) != 1)
1463  std::cout << "Error 1 in row " << i << std::endl;
1464  if (P(i,j) == 1 && set)
1465  std::cout << "Error 2 in row " << i << std::endl;
1466  if (P(i,j) == 1 && !set)
1467  set = true;
1468  }
1469  }
1470  }
1471 
1472  if (Pointvector.is_fpoint(i))
1473  for (unsigned int j=0; j<P.size2(); ++j)
1474  {
1475  if (P.isnonzero(i,j) && j> Pointvector.get_cpoints()-1)
1476  std::cout << "Error 3 in row " << i << std::endl;
1477  if (P.isnonzero(i,j))
1478  {
1479  bool set = false;
1480  for (unsigned int k=0; k<P.size1(); ++k)
1481  {
1482  if (P.isnonzero(k,j))
1483  {
1484  if (Pointvector.is_cpoint(k) && P(k,j) == 1 && A.isnonzero(i,k))
1485  set = true;
1486  }
1487  }
1488  if (!set)
1489  std::cout << "Error 4 in row " << i << std::endl;
1490  }
1491  }
1492  }
1493  }
1494 
1495 
1496  } //namespace amg
1497  }
1498  }
1499 }
1500 
1501 #endif
void clear_influencelists()
Definition: amg_base.hpp:991
iterator begin()
Definition: amg_base.hpp:352
void resize(vcl_size_t i, vcl_size_t j, bool preserve=true)
Definition: adapter.hpp:377
ConstAdapterType::const_iterator2 const_iterator2
Definition: amg_base.hpp:409
amg_sparsevector_iterator(InternalType &vec, bool begin=true)
The constructor.
Definition: amg_base.hpp:218
Debug functionality for AMG. To be removed.
std::size_t vcl_size_t
Definition: forwards.h:58
double get_jacobiweight() const
Definition: amg_base.hpp:102
void set_offset(unsigned int offset)
Definition: amg_base.hpp:850
void set_coarselevels(int coarselevels)
Definition: amg_base.hpp:110
void switch_ftoc()
Definition: amg_base.hpp:904
iterator end()
Definition: amg_base.hpp:354
void set_as(double jacobiweight)
Definition: amg_base.hpp:98
ScalarType value_type
Definition: amg_base.hpp:257
const_iterator2 end2(bool trans=false) const
Definition: amg_base.hpp:768
void add_influencing_point(amg_point *point)
Definition: amg_base.hpp:867
const_iterator2 begin2() const
Definition: adapter.hpp:201
vcl_size_t size1()
Definition: amg_base.hpp:653
amg_sparsematrix(MatrixType const &mat)
Constructor. Builds matrix via another matrix type. (Only necessary feature of this other matrix type...
Definition: amg_base.hpp:459
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:330
void addscalar(IteratorType &iter, unsigned int i, unsigned int, ScalarType s)
Definition: amg_base.hpp:312
A class for the AMG points. Saves point index and influence measure Holds information whether point i...
Definition: amg_base.hpp:816
std::vector< std::map< unsigned int, ScalarType > > * get_internal_pointer()
Definition: amg_base.hpp:776
void addscalar(IteratorType &iter, unsigned int i, unsigned int j, ScalarType s)
Definition: amg_base.hpp:566
amg_nonzero_scalar(InternalType *m, IteratorType &iter, unsigned int i, unsigned int j, ScalarType s=0)
The constructor.
Definition: amg_base.hpp:142
unsigned int get_influence() const
Definition: amg_base.hpp:854
ScalarType operator=(const ScalarType value)
Assignment operator. Writes value into matrix at the given position.
Definition: amg_base.hpp:151
void set_coarse_index(unsigned int index)
Definition: amg_base.hpp:878
void update_cf(amg_point *point)
Definition: amg_base.hpp:982
bool get_trans() const
Definition: amg_base.hpp:520
void get_F(VectorType &vec) const
Definition: amg_base.hpp:1139
const_iterator begin() const
Definition: amg_base.hpp:1003
void get_sorting(VectorType &vec) const
Definition: amg_base.hpp:1114
amg_point * get_nextpoint()
Definition: amg_base.hpp:1025
void test_triplematprod(SparseMatrixType &A, SparseMatrixType &P, SparseMatrixType &A_i1)
Test triple-matrix product by comparing it to ublas functions. Very slow for large matrices! ...
Definition: amg_base.hpp:1417
unsigned int levels_
Definition: amg_base.hpp:1181
iterator begin()
Definition: amg_base.hpp:1001
boost::numeric::ublas::vector< boost::numeric::ublas::vector< unsigned int > > Offset
Definition: amg_base.hpp:1178
A class for the sparse vector type.
Definition: amg_base.hpp:254
void clear()
Definition: amg_base.hpp:644
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 calc_influence()
Definition: amg_base.hpp:881
void init(unsigned int levels, unsigned int threads=0)
Definition: amg_base.hpp:1183
void slice(unsigned int level, InternalType1 const &A, InternalType2 const &Pointvector)
Definition: amg_base.hpp:1213
unsigned int get_coarse() const
Definition: amg_base.hpp:90
iterator2 end2(bool trans=false)
Definition: amg_base.hpp:731
const_iterator1 end1() const
Definition: adapter.hpp:196
void resize(unsigned int i, unsigned int j, bool preserve=true)
Definition: amg_base.hpp:634
void add(unsigned int i, ScalarType s)
Definition: amg_base.hpp:293
void clear_influencing()
Definition: amg_base.hpp:872
iterator1 begin1()
Definition: adapter.hpp:363
unsigned int get_coarse_index() const
Definition: amg_base.hpp:877
unsigned int get_postsmooth() const
Definition: amg_base.hpp:108
bool is_undecided() const
Definition: amg_base.hpp:860
void printmatrix(MatrixType &, int)
Definition: amg_debug.hpp:77
iterator2 begin2(bool trans=false)
Definition: amg_base.hpp:716
void set_aggregate(unsigned int aggregate)
Definition: amg_base.hpp:855
void set_index(unsigned int index)
Definition: amg_base.hpp:852
VectorType::iterator iterator
Definition: amg_base.hpp:948
void add_influenced_point(amg_point *point)
Definition: amg_base.hpp:869
void clear_cf()
Definition: amg_base.hpp:988
amg_point * operator[](unsigned int i) const
Definition: amg_base.hpp:1000
unsigned int get_cpoints() const
Definition: amg_base.hpp:1014
unsigned int index() const
Definition: amg_base.hpp:247
NonzeroScalarType operator[](unsigned int i)
Definition: amg_base.hpp:330
unsigned int get_offset()
Definition: amg_base.hpp:851
bool operator==(self_type other)
Definition: amg_base.hpp:226
amg_sparsematrix()
Standard constructor.
Definition: amg_base.hpp:412
unsigned int internal_size() const
Definition: amg_base.hpp:286
const_iterator begin_influenced() const
Definition: amg_base.hpp:913
AdapterType::iterator1 iterator1
Definition: amg_base.hpp:406
amg_point(unsigned int index, unsigned int size)
The constructor.
Definition: amg_base.hpp:844
iterator1 end1()
Definition: adapter.hpp:364
void make_cpoint()
Definition: amg_base.hpp:890
A const iterator for sparse matrices of type std::vector > ...
Definition: adapter.hpp:48
void set_interpolweight(double interpolweight)
Definition: amg_base.hpp:101
ListType::const_iterator const_iterator
Definition: amg_base.hpp:840
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
const_iterator end_influencing() const
Definition: amg_base.hpp:910
unsigned int get_aggregate()
Definition: amg_base.hpp:856
void init_points()
Definition: amg_base.hpp:961
unsigned int number_influencing() const
Definition: amg_base.hpp:863
void sort()
Definition: amg_base.hpp:1019
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:91
void resize(unsigned int size)
Definition: amg_base.hpp:1006
ListType::iterator iterator
Definition: amg_base.hpp:839
void add(unsigned int i, unsigned int j, ScalarType s)
Definition: amg_base.hpp:542
void set_presmooth(int presmooth)
Definition: amg_base.hpp:104
self_type & operator--() const
Definition: amg_base.hpp:243
boost::numeric::ublas::vector< InternalType1 > A_slice
Definition: amg_base.hpp:1175
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 make_cpoint(amg_point *point)
Definition: amg_base.hpp:1056
boost::numeric::ublas::vector< InternalType2 > Pointvector_slice
Definition: amg_base.hpp:1176
const_iterator end() const
Definition: amg_base.hpp:355
bool operator!=(self_type other)
Definition: amg_base.hpp:233
void make_fpoint(amg_point *point)
Definition: amg_base.hpp:1063
void clear_influenced()
Definition: amg_base.hpp:874
iterator begin_influencing()
Definition: amg_base.hpp:907
unsigned int get_fpoints() const
Definition: amg_base.hpp:1016
vcl_size_t size2()
Definition: amg_base.hpp:670
ScalarType operator++()
Definition: amg_base.hpp:189
void test_interpolation(SparseMatrixType &A, SparseMatrixType &P, PointVectorType &Pointvector)
Test if interpolation matrix makes sense. Only vanilla test though! Only checks if basic requirements...
Definition: amg_base.hpp:1451
void removescalar(IteratorType &iter, unsigned int)
Definition: amg_base.hpp:327
bool is_cpoint() const
Definition: amg_base.hpp:858
void get_influence(VectorType &vec) const
Definition: amg_base.hpp:1105
iterator begin_influenced()
Definition: amg_base.hpp:911
void do_trans()
Definition: amg_base.hpp:488
void clear()
Definition: amg_base.hpp:288
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 add_influence(amg_point *point, unsigned int add)
Definition: amg_base.hpp:1038
unsigned int get_coarselevels() const
Definition: amg_base.hpp:111
void set_threshold(double threshold)
Definition: amg_base.hpp:95
amg_sparsematrix(unsigned int i, unsigned int j)
Constructor. Builds matrix of size (i,j).
Definition: amg_base.hpp:422
Adapts a constant sparse matrix type made up from std::vector > to bas...
Definition: adapter.hpp:176
A non-const iterator for sparse matrices of type std::vector > ...
Definition: adapter.hpp:230
const_iterator begin_influencing() const
Definition: amg_base.hpp:909
iterator1 end1(bool trans=false)
Definition: amg_base.hpp:701
void resize(unsigned int size)
Definition: amg_base.hpp:282
void make_fpoint()
Definition: amg_base.hpp:897
AdapterType::iterator2 iterator2
Definition: amg_base.hpp:407
bool isnonzero(unsigned int i, unsigned int j) const
Definition: amg_base.hpp:523
A class for the sparse matrix type. Uses vector of maps as data structure for higher performance and ...
Definition: amg_base.hpp:376
iterator1 begin1(bool trans=false)
Definition: amg_base.hpp:686
unsigned int size() const
Definition: amg_base.hpp:283
ScalarType operator()(unsigned int i, unsigned int j) const
Definition: amg_base.hpp:612
void set_postsmooth(int postsmooth)
Definition: amg_base.hpp:107
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_threshold() const
Definition: amg_base.hpp:96
ScalarType & operator*() const
Definition: amg_base.hpp:245
double get_interpolweight() const
Definition: amg_base.hpp:99
const_iterator1 begin1() const
Definition: amg_base.hpp:746
unsigned int threads_
Definition: amg_base.hpp:1180
void get_Aggregates(MatrixType &mat) const
Definition: amg_base.hpp:1151
NonzeroScalarType operator()(unsigned int i, unsigned int j)
Definition: amg_base.hpp:589
unsigned int get_index() const
Definition: amg_base.hpp:853
void set_interpol(unsigned int interpol)
Definition: amg_base.hpp:92
self_type & operator++() const
Definition: amg_base.hpp:241
amg_sparsevector_iterator< InternalType > iterator
Definition: amg_base.hpp:270
A class for a scalar that can be written to the sparse matrix or sparse vector datatypes.
Definition: amg_base.hpp:124
void build_index()
Definition: amg_base.hpp:1078
void clear()
Definition: adapter.hpp:388
InternalType::const_iterator const_iterator
Definition: amg_base.hpp:271
unsigned int get_interpol() const
Definition: amg_base.hpp:93
void get_influence_matrix(MatrixType &mat) const
Definition: amg_base.hpp:1095
void amg_mat_prod(SparseMatrixType &A, SparseMatrixType &B, SparseMatrixType &RES)
Sparse matrix product. Calculates RES = A*B.
Definition: amg_base.hpp:1319
void delete_points()
Definition: amg_base.hpp:967
const_iterator end_influenced() const
Definition: amg_base.hpp:914
const_iterator end() const
Definition: amg_base.hpp:1004
ConstAdapterType::const_iterator1 const_iterator1
Definition: amg_base.hpp:408
vcl_size_t size1() const
Definition: amg_base.hpp:661
ScalarType operator+=(const ScalarType value)
Addition operator. Adds a constant.
Definition: amg_base.hpp:164
iterator end()
Definition: amg_base.hpp:1002
bool operator()(amg_point *l, amg_point *r) const
Definition: amg_base.hpp:922
bool is_influencing(amg_point *point) const
Definition: amg_base.hpp:865
void set_coarse(unsigned int coarse)
Definition: amg_base.hpp:89
void switch_ftoc(amg_point *point)
Definition: amg_base.hpp:1070
void join(unsigned int level, InternalType2 &Pointvector) const
Definition: amg_base.hpp:1226
const_iterator begin() const
Definition: amg_base.hpp:353
amg_sparsematrix(std::vector< std::map< unsigned int, ScalarType > > const &mat)
Constructor. Builds matrix via std::vector by copying memory (Only necessary feature of thi...
Definition: amg_base.hpp:440
Defines an iterator for the sparse vector type.
Definition: amg_base.hpp:203
iterator2 begin2()
Definition: adapter.hpp:369
void removescalar(IteratorType &iter, unsigned int i)
Definition: amg_base.hpp:582
const_iterator1 begin1() const
Definition: adapter.hpp:195
Adapts a non-const sparse matrix type made up from std::vector > to ba...
Definition: adapter.hpp:345
const_iterator2 begin2(bool trans=false) const
Definition: amg_base.hpp:761
const_iterator1 end1(bool trans=false) const
Definition: amg_base.hpp:754
iterator end_influenced()
Definition: amg_base.hpp:912
amg_pointvector(unsigned int size=0)
The constructor.
Definition: amg_base.hpp:954
unsigned int add_influence(int add)
Definition: amg_base.hpp:884
unsigned int index()
Definition: amg_base.hpp:248
bool is_fpoint() const
Definition: amg_base.hpp:859
const_iterator2 end2() const
Definition: adapter.hpp:202
Comparison class for the sorted set of points in amg_pointvector. Set is sorted by influence measure ...
Definition: amg_base.hpp:919
iterator end_influencing()
Definition: amg_base.hpp:908
void get_C(VectorType &vec) const
Definition: amg_base.hpp:1127
void add_point(amg_point *point)
Definition: amg_base.hpp:973
void set_trans(bool mode)
Definition: amg_base.hpp:513
bool isnonzero(unsigned int i) const
Definition: amg_base.hpp:358
unsigned int size() const
Definition: amg_base.hpp:1011
amg_tag(unsigned int coarse=1, unsigned int interpol=1, double threshold=0.25, double interpolweight=0.2, double jacobiweight=1, unsigned int presmooth=1, unsigned int postsmooth=1, unsigned int coarselevels=0)
The constructor.
Definition: amg_base.hpp:76
iterator2 end2()
Definition: adapter.hpp:370
amg_sparsevector(unsigned int size=0)
The constructor.
Definition: amg_base.hpp:277
A class for the matrix slicing for parallel coarsening schemes (RS0/RS3).
Definition: amg_base.hpp:1168
VectorType::const_iterator const_iterator
Definition: amg_base.hpp:949
vcl_size_t size2() const
Definition: amg_base.hpp:678
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
ScalarType value_type
Definition: amg_base.hpp:379