ViennaCL - The Vienna Computing Library  1.2.0
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-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 
26 #include <boost/numeric/ublas/vector.hpp>
27 #include <cmath>
28 #include <set>
29 #include <list>
30 #include <algorithm>
31 
32 #include <map>
33 #ifdef _OPENMP
34 #include <omp.h>
35 #endif
36 
37 #include "amg_debug.hpp"
38 
39 #define VIENNACL_AMG_COARSE_RS 1
40 #define VIENNACL_AMG_COARSE_ONEPASS 2
41 #define VIENNACL_AMG_COARSE_RS0 3
42 #define VIENNACL_AMG_COARSE_RS3 4
43 #define VIENNACL_AMG_COARSE_AG 5
44 #define VIENNACL_AMG_INTERPOL_DIRECT 1
45 #define VIENNACL_AMG_INTERPOL_CLASSIC 2
46 #define VIENNACL_AMG_INTERPOL_AG 3
47 #define VIENNACL_AMG_INTERPOL_SA 4
48 
49 namespace viennacl
50 {
51  namespace linalg
52  {
53  namespace detail
54  {
55  namespace amg
56  {
59  class amg_tag
60  {
61  public:
74  amg_tag(unsigned int coarse = 1,
75  unsigned int interpol = 1,
76  double threshold = 0.25,
77  double interpolweight = 0.2,
78  double jacobiweight = 1,
79  unsigned int presmooth = 1,
80  unsigned int postsmooth = 1,
81  unsigned int coarselevels = 0)
82  : _coarse(coarse), _interpol(interpol),
83  _threshold(threshold), _interpolweight(interpolweight), _jacobiweight(jacobiweight),
84  _presmooth(presmooth), _postsmooth(postsmooth), _coarselevels(coarselevels) {};
85 
86  // Getter-/Setter-Functions
87  void set_coarse(unsigned int coarse) { if (coarse > 0) _coarse = coarse; }
88  unsigned int get_coarse() const { return _coarse; }
89 
90  void set_interpol(unsigned int interpol) { if (interpol > 0) _interpol = interpol; }
91  unsigned int get_interpol() const { return _interpol; }
92 
93  void set_threshold(double threshold) { if (threshold > 0 && threshold <= 1) _threshold = threshold; }
94  double get_threshold() const{ return _threshold; }
95 
96  void set_as(double jacobiweight) { if (jacobiweight > 0 && jacobiweight <= 2) _jacobiweight = jacobiweight; }
97  double get_interpolweight() const { return _interpolweight; }
98 
99  void set_interpolweight(double interpolweight) { if (interpolweight > 0 && interpolweight <= 2) _interpolweight = interpolweight; }
100  double get_jacobiweight() const { return _jacobiweight; }
101 
102  void set_presmooth(int presmooth) { if (presmooth >= 0) _presmooth = presmooth; }
103  unsigned int get_presmooth() const { return _presmooth; }
104 
105  void set_postsmooth(int postsmooth) { if (postsmooth >= 0) _postsmooth = postsmooth; }
106  unsigned int get_postsmooth() const { return _postsmooth; }
107 
108  void set_coarselevels(int coarselevels) { if (coarselevels >= 0) _coarselevels = coarselevels; }
109  unsigned int get_coarselevels() const { return _coarselevels; }
110 
111  private:
112  unsigned int _coarse, _interpol;
113  double _threshold, _interpolweight, _jacobiweight;
114  unsigned int _presmooth, _postsmooth, _coarselevels;
115  };
116 
121  template <typename InternalType, typename IteratorType, typename ScalarType>
123  {
124  private:
125  InternalType *_m;
126  IteratorType _iter;
127  unsigned int _i,_j;
128  ScalarType _s;
129 
130  public:
132 
140  amg_nonzero_scalar(InternalType *m,
141  IteratorType & iter,
142  unsigned int i,
143  unsigned int j,
144  ScalarType s = 0): _m(m), _iter(iter), _i(i), _j(j), _s(s) {}
145 
149  ScalarType operator = (const ScalarType value)
150  {
151  _s = value;
152  // Only write if scalar is nonzero
153  if (_s == 0) return _s;
154  // Write to _m using iterator _iter or indices (_i,_j)
155  _m->addscalar (_iter,_i,_j,_s);
156  return _s;
157  }
158 
162  ScalarType operator += (const ScalarType value)
163  {
164  // If zero is added, then no change necessary
165  if (value == 0)
166  return _s;
167 
168  _s += value;
169  // Remove entry if resulting scalar is zero
170  if (_s == 0)
171  {
172  _m->removescalar(_iter,_i);
173  return _s;
174  }
175  //Write to _m using iterator _iter or indices (_i,_j)
176  _m->addscalar (_iter,_i,_j,_s);
177  return _s;
178  }
179  ScalarType operator ++ (int)
180  {
181  _s++;
182  if (_s == 0)
183  _m->removescalar(_iter,_i);
184  _m->addscalar (_iter,_i,_j,_s);
185  return _s;
186  }
187  ScalarType operator ++ ()
188  {
189  _s++;
190  if (_s == 0)
191  _m->removescalar(_iter,_i);
192  _m->addscalar (_iter,_i,_j,_s);
193  return _s;
194  }
195  operator ScalarType (void) { return _s; }
196  };
197 
200  template <typename InternalType>
202  {
203  private:
205  typedef typename InternalType::mapped_type ScalarType;
206 
207  InternalType & internal_vec;
208  typename InternalType::iterator iter;
209 
210  public:
211 
216  amg_sparsevector_iterator(InternalType & vec, bool begin=true): internal_vec(vec)
217  {
218  if (begin)
219  iter = internal_vec.begin();
220  else
221  iter = internal_vec.end();
222  }
223 
225  {
226  if (iter == other.iter)
227  return true;
228  else
229  return false;
230  }
232  {
233  if (iter != other.iter)
234  return true;
235  else
236  return false;
237  }
238 
239  self_type & operator ++ () const { iter++; return *this; }
240  self_type & operator ++ () { iter++; return *this; }
241  self_type & operator -- () const { iter--; return *this; }
242  self_type & operator -- () { iter--; return *this; }
243  ScalarType & operator * () const { return (*iter).second; }
244  ScalarType & operator * () { return (*iter).second; }
245  unsigned int index() const { return (*iter).first; }
246  unsigned int index() { return (*iter).first; }
247  };
248 
251  template <typename ScalarType>
253  {
254  public:
255  typedef ScalarType value_type;
256 
257  private:
258  // A map is used internally which saves all non-zero elements with pairs of (index,value)
259  typedef std::map<unsigned int,ScalarType> InternalType;
262 
263  // Size is only a dummy variable. Not needed for internal map structure but for compatible vector interface.
264  unsigned int _size;
265  InternalType internal_vector;
266 
267  public:
269  typedef typename InternalType::const_iterator const_iterator;
270 
271  public:
275  amg_sparsevector(unsigned int size = 0): _size(size)
276  {
277  internal_vector = InternalType();
278  }
279 
280  void resize(unsigned int size) { _size = size; }
281  unsigned int size() const { return _size;}
282 
283  // Returns number of non-zero entries in vector equal to the size of the underlying map.
284  unsigned int internal_size() const { return internal_vector.size(); }
285  // Delete underlying map.
286  void clear() { internal_vector.clear(); }
287  // Remove entry at position i.
288  void remove(unsigned int i) { internal_vector.erase(i); }
289 
290  // Add s to the entry at position i
291  void add (unsigned int i, ScalarType s)
292  {
293  typename InternalType::iterator iter = internal_vector.find(i);
294  // If there is no entry at position i, add new entry at that position
295  if (iter == internal_vector.end())
296  addscalar(iter,i,i,s);
297  else
298  {
299  s += (*iter).second;
300  // If new value is zero, then erase the entry, otherwise write new value
301  if (s != 0)
302  (*iter).second = s;
303  else
304  internal_vector.erase(iter);
305  }
306  }
307 
308  // Write to the map. Is called from non-zero scalar type.
309  template <typename IteratorType>
310  void addscalar(IteratorType & iter, unsigned int i, unsigned int j, ScalarType s)
311  {
312  // Don't write if value is zero
313  if (s == 0)
314  return;
315 
316  // If entry is already present, overwrite value, otherwise make new entry
317  if (iter != internal_vector.end())
318  (*iter).second = s;
319  else
320  internal_vector[i] = s;
321  }
322 
323  // Remove value from the map. Is called from non-zero scalar type.
324  template <typename IteratorType>
325  void removescalar(IteratorType & iter, unsigned int i) { internal_vector.erase(iter); }
326 
327  // Bracket operator. Returns non-zero scalar type with actual values of the respective entry which calls addscalar/removescalar after value is altered.
329  {
330  typename InternalType::iterator it = internal_vector.find(i);
331  // If value is already present then build non-zero scalar with actual value, otherwise 0.
332  if (it != internal_vector.end())
333  return NonzeroScalarType (this,it,i,i,(*it).second);
334  else
335  return NonzeroScalarType (this,it,i,i,0);
336  }
337 
338  // Use internal data structure directly for read-only access. No need to use non-zero scalar as no write access possible.
339  ScalarType operator [] (unsigned int i) const
340  {
341  const_iterator it = internal_vector.find(i);
342 
343  if (it != internal_vector.end())
344  return (*it).second;
345  else
346  return 0;
347  }
348 
349  // Iterator functions.
350  iterator begin() { return iterator(internal_vector); }
351  const_iterator begin() const { return internal_vector.begin(); }
352  iterator end() { return iterator(internal_vector,false); }
353  const_iterator end() const { return internal_vector.end(); }
354 
355  // checks whether value at index i is nonzero. More efficient than doing [] == 0.
356  bool isnonzero(unsigned int i) const { return internal_vector.find(i) != internal_vector.end(); }
357 
358  // Copies data into a ublas vector type.
359  operator boost::numeric::ublas::vector<ScalarType> (void)
360  {
361  boost::numeric::ublas::vector<ScalarType> vec (_size);
362  for (iterator iter = begin(); iter != end(); ++iter)
363  vec [iter.index()] = *iter;
364  return vec;
365  }
366  };
367 
373  template <typename ScalarType>
375  {
376  public:
377  typedef ScalarType value_type;
378  private:
379  typedef std::map<unsigned int,ScalarType> RowType;
380  typedef std::vector<RowType> InternalType;
382 
383  // Adapter is used for certain functionality, especially iterators.
386 
387  // Non-zero scalar is used to write to the matrix.
389 
390  // Holds matrix coefficients.
391  InternalType internal_mat;
392  // Holds matrix coefficient of transposed matrix if built.
393  // Note: Only internal_mat is written using operators and methods while internal_mat_trans is built from internal_mat using do_trans().
394  InternalType internal_mat_trans;
395  // Saves sizes.
396  size_t s1, s2;
397 
398  // True if the transposed of the matrix is used (for calculations, iteration, etc.).
399  bool transposed_mode;
400  // True if the transposed is already built (saved in internal_mat_trans) and also up to date (no changes to internal_mat).
401  bool transposed;
402 
403  public:
408 
411  {
412  transposed_mode = false;
413  transposed = false;
414  }
415 
420  amg_sparsematrix (unsigned int i, unsigned int j)
421  {
422  AdapterType a (internal_mat);
423  a.resize (i,j,false);
424  AdapterType a_trans (internal_mat_trans);
425  a_trans.resize (j,i,false);
426  s1 = i;
427  s2 = j;
428  a.clear();
429  a_trans.clear();
430  transposed_mode = false;
431  transposed = false;
432  }
433 
438  amg_sparsematrix (std::vector<std::map<unsigned int, ScalarType> > const & mat)
439  {
440  AdapterType a (internal_mat);
441  AdapterType a_trans (internal_mat_trans);
442  a.resize(mat.size(), mat.size());
443  a_trans.resize(mat.size(), mat.size());
444 
445  internal_mat = mat;
446  s1 = s2 = mat.size();
447 
448  transposed_mode = false;
449  transposed = false;
450  }
451 
456  template <typename MatrixType>
457  amg_sparsematrix (MatrixType const & mat)
458  {
459  AdapterType a (internal_mat);
460  AdapterType a_trans (internal_mat_trans);
461  a.resize(mat.size1(), mat.size2());
462  a_trans.resize (mat.size2(), mat.size1());
463  s1 = mat.size1();
464  s2 = mat.size2();
465  a.clear();
466  a_trans.clear();
467 
468  for (typename MatrixType::const_iterator1 row_iter = mat.begin1(); row_iter != mat.end1(); ++row_iter)
469  {
470  for (typename MatrixType::const_iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
471  {
472  if (*col_iter != 0)
473  {
474  unsigned int x = col_iter.index1();
475  unsigned int y = col_iter.index2();
476  a (x,y) = *col_iter;
477  a_trans (y,x) = *col_iter;
478  }
479  }
480  }
481  transposed_mode = false;
482  transposed = true;
483  }
484 
485  // Build transposed of the current matrix.
486  void do_trans()
487  {
488  // Do it only once if called in a parallel section
489  #ifdef _OPENMP
490  #pragma omp critical
491  #endif
492  {
493  // Only build transposed if it is not built or not up to date
494  if (!transposed)
495  {
496  // Mode has to be set to standard mode temporarily
497  bool save_mode = transposed_mode;
498  transposed_mode = false;
499 
500  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
501  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
502  internal_mat_trans[col_iter.index2()][col_iter.index1()] = *col_iter;
503 
504  transposed_mode = save_mode;
505  transposed = true;
506  }
507  }
508  } //do_trans()
509 
510  // Set transposed mode (true=transposed, false=regular)
511  void set_trans(bool mode)
512  {
513  transposed_mode = mode;
514  if (mode)
515  do_trans();
516  }
517 
518  bool get_trans() const { return transposed_mode; }
519 
520  // Checks whether coefficient (i,j) is non-zero. More efficient than using (i,j) == 0.
521  bool isnonzero (unsigned int i, unsigned int j) const
522  {
523  if (!transposed_mode)
524  {
525  if (internal_mat[i].find(j) != internal_mat[i].end())
526  return true;
527  else
528  return false;
529  }
530  else
531  {
532  if (internal_mat[j].find(i) != internal_mat[j].end())
533  return true;
534  else
535  return false;
536  }
537  } //isnonzero()
538 
539  // Add s to value at (i,j)
540  void add (unsigned int i, unsigned int j, ScalarType s)
541  {
542  // If zero is added then do nothing.
543  if (s == 0)
544  return;
545 
546  typename RowType::iterator col_iter = internal_mat[i].find(j);
547  // If there is no entry at position (i,j), then make new entry.
548  if (col_iter == internal_mat[i].end())
549  addscalar(col_iter,i,j,s);
550  else
551  {
552  s += (*col_iter).second;
553  // Update value and erase entry if value is zero.
554  if (s != 0)
555  (*col_iter).second = s;
556  else
557  internal_mat[i].erase(col_iter);
558  }
559  transposed = false;
560  } //add()
561 
562  // Write to the internal data structure. Is called from non-zero scalar type.
563  template <typename IteratorType>
564  void addscalar(IteratorType & iter, unsigned int i, unsigned int j, ScalarType s)
565  {
566  // Don't write if value is zero
567  if (s == 0)
568  return;
569 
570  if (iter != internal_mat[i].end())
571  (*iter).second = s;
572  else
573  internal_mat[i][j] = s;
574 
575  transposed = false;
576  }
577 
578  // Remove entry from internal data structure. Is called from non-zero scalar type.
579  template <typename IteratorType>
580  void removescalar(IteratorType & iter, unsigned int i)
581  {
582  internal_mat[i].erase(iter);
583  transposed = false;
584  }
585 
586  // Return non-zero scalar at position (i,j). Value is written to the non-zero scalar and updated via addscalar()/removescalar().
587  NonzeroScalarType operator()(unsigned int i, unsigned int j)
588  {
589  typename RowType::iterator iter;
590 
591  if (!transposed_mode)
592  {
593  iter = internal_mat[i].find(j);
594  if (iter != internal_mat[i].end())
595  return NonzeroScalarType (this,iter,i,j,(*iter).second);
596  else
597  return NonzeroScalarType (this,iter,i,j,0);
598  }
599  else
600  {
601  iter = internal_mat[j].find(i);
602  if (iter != internal_mat[j].end())
603  return NonzeroScalarType (this,iter,j,i,(*iter).second);
604  else
605  return NonzeroScalarType (this,iter,j,i,0);
606  }
607  }
608 
609  // For read-only access return the actual value directly. Non-zero datatype not needed as no write access possible.
610  ScalarType operator()(unsigned int i, unsigned int j) const
611  {
612  typename RowType::const_iterator iter;
613 
614  if (!transposed_mode)
615  {
616  iter = internal_mat[i].find(j);
617  if (iter != internal_mat[i].end())
618  return (*iter).second;
619  else
620  return 0;
621  }
622  else
623  {
624  iter = internal_mat[j].find(i);
625  if (iter != internal_mat[j].end())
626  return (*iter).second;
627  else
628  return 0;
629  }
630  }
631 
632  void resize(unsigned int i, unsigned int j, bool preserve = true)
633  {
634  AdapterType a (internal_mat);
635  a.resize(i,j,preserve);
636  AdapterType a_trans (internal_mat_trans);
637  a_trans.resize(j,i,preserve);
638  s1 = i;
639  s2 = j;
640  }
641 
642  void clear()
643  {
644  AdapterType a (internal_mat);
645  a.clear();
646  AdapterType a_trans (internal_mat_trans);
647  a_trans.clear();
648  transposed = true;
649  }
650 
651  size_t size1()
652  {
653  if (!transposed_mode)
654  return s1;
655  else
656  return s2;
657  }
658 
659  size_t size1() const
660  {
661  if (!transposed_mode)
662  return s1;
663  else
664  return s2;
665  }
666 
667 
668  size_t size2()
669  {
670  if (!transposed_mode)
671  return s2;
672  else
673  return s1;
674  }
675  size_t size2() const
676  {
677  if (!transposed_mode)
678  return s2;
679  else
680  return s1;
681  }
682 
683  iterator1 begin1(bool trans = false)
684  {
685  if (!trans && !transposed_mode)
686  {
687  AdapterType a (internal_mat);
688  return a.begin1();
689  }
690  else
691  {
692  do_trans();
693  AdapterType a_trans (internal_mat_trans);
694  return a_trans.begin1();
695  }
696  }
697 
698  iterator1 end1(bool trans = false)
699  {
700  if (!trans && !transposed_mode)
701  {
702  AdapterType a (internal_mat);
703  return a.end1();
704  }
705  else
706  {
707  //do_trans();
708  AdapterType a_trans (internal_mat_trans);
709  return a_trans.end1();
710  }
711  }
712 
713  iterator2 begin2(bool trans = false)
714  {
715  if (!trans && !transposed_mode)
716  {
717  AdapterType a (internal_mat);
718  return a.begin2();
719  }
720  else
721  {
722  do_trans();
723  AdapterType a_trans (internal_mat_trans);
724  return a_trans.begin2();
725  }
726  }
727 
728  iterator2 end2(bool trans = false)
729  {
730  if (!trans && !transposed_mode)
731  {
732  AdapterType a (internal_mat);
733  return a.end2();
734  }
735  else
736  {
737  //do_trans();
738  AdapterType a_trans (internal_mat_trans);
739  return a_trans.end2();
740  }
741  }
742 
744  {
745  // Const_iterator of transposed can only be used if transposed matrix is already built and up to date.
746  assert((!transposed_mode || (transposed_mode && transposed)) && "Error: Cannot build const_iterator when transposed has not been built yet!");
747  ConstAdapterType a_const (internal_mat);
748  return a_const.begin1();
749  }
750 
751  const_iterator1 end1(bool trans = false) const
752  {
753  assert((!transposed_mode || (transposed_mode && transposed)) && "Error: Cannot build const_iterator when transposed has not been built yet!");
754  ConstAdapterType a_const (internal_mat);
755  return a_const.end1();
756  }
757 
758  const_iterator2 begin2(bool trans = false) const
759  {
760  assert((!transposed_mode || (transposed_mode && transposed)) && "Error: Cannot build const_iterator when transposed has not been built yet!");
761  ConstAdapterType a_const (internal_mat);
762  return a_const.begin2();
763  }
764 
765  const_iterator2 end2(bool trans = false) const
766  {
767  assert((!transposed_mode || (transposed_mode && transposed)) && "Error: Cannot build const_iterator when transposed has not been built yet!");
768  ConstAdapterType a_const (internal_mat);
769  return a_const.end2();
770  }
771 
772  // Returns pointer to the internal data structure. Improves performance of copy operation to GPU.
773  std::vector<std::map<unsigned int, ScalarType> > * get_internal_pointer()
774  {
775  if (!transposed_mode)
776  return &internal_mat;
777 
778  if (!transposed)
779  do_trans();
780  return &internal_mat_trans;
781  }
782  operator boost::numeric::ublas::compressed_matrix<ScalarType> (void)
783  {
784  boost::numeric::ublas::compressed_matrix<ScalarType> mat;
785  mat.resize(size1(),size2(),false);
786  mat.clear();
787 
788  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
789  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
790  mat (col_iter.index1(), col_iter.index2()) = *col_iter;
791 
792  return mat;
793  }
794  operator boost::numeric::ublas::matrix<ScalarType> (void)
795  {
796  boost::numeric::ublas::matrix<ScalarType> mat;
797  mat.resize(size1(),size2(),false);
798  mat.clear();
799 
800  for (iterator1 row_iter = begin1(); row_iter != end1(); ++row_iter)
801  for (iterator2 col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
802  mat (col_iter.index1(), col_iter.index2()) = *col_iter;
803 
804  return mat;
805  }
806  };
807 
813  class amg_point
814  {
815  private:
817 
818  unsigned int _index;
819  unsigned int _influence;
820  // Determines whether point is undecided.
821  bool _undecided;
822  // Determines wheter point is C point (true) or F point (false).
823  bool _cpoint;
824  unsigned int _coarse_index;
825  // 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
826  unsigned int _offset;
827  // Aggregate the point belongs to.
828  unsigned int _aggregate;
829 
830  // Holds all points influencing this point.
831  ListType influencing_points;
832  // Holds all points that are influenced by this point.
833  ListType influenced_points;
834 
835  public:
838 
841  amg_point (unsigned int index, unsigned int size): _index(index), _influence(0), _undecided(true), _cpoint(false), _coarse_index(0), _offset(0), _aggregate(0)
842  {
843  influencing_points = ListType(size);
844  influenced_points = ListType(size);
845  }
846 
847  void set_offset(unsigned int offset) { _offset = offset; }
848  unsigned int get_offset() { return _offset; }
849  void set_index(unsigned int index) { _index = index+_offset; }
850  unsigned int get_index() const { return _index-_offset; }
851  unsigned int get_influence() const { return _influence; }
852  void set_aggregate(unsigned int aggregate) { _aggregate = aggregate; }
853  unsigned int get_aggregate () { return _aggregate; }
854 
855  bool is_cpoint() const { return _cpoint && !_undecided; }
856  bool is_fpoint() const { return !_cpoint && !_undecided; }
857  bool is_undecided() const { return _undecided; }
858 
859  // Returns number of influencing points
860  unsigned int number_influencing() const { return influencing_points.internal_size(); }
861  // Returns true if *point is influencing this point
862  bool is_influencing(amg_point* point) const { return influencing_points.isnonzero(point->get_index()+_offset); }
863  // Add *point to influencing points
864  void add_influencing_point(amg_point* point) { influencing_points[point->get_index()+_offset] = point; }
865  // Add *point to influenced points
866  void add_influenced_point(amg_point* point) { influenced_points[point->get_index()+_offset] = point; }
867 
868  // Clear influencing points
869  void clear_influencing() { influencing_points.clear(); }
870  // Clear influenced points
871  void clear_influenced() {influenced_points.clear(); }
872 
873 
874  unsigned int get_coarse_index() const { return _coarse_index; }
875  void set_coarse_index(unsigned int index) { _coarse_index = index; }
876 
877  // Calculates the initial influence measure equal to the number of influenced points.
878  void calc_influence() { _influence = influenced_points.internal_size(); }
879 
880  // Add to influence measure.
881  unsigned int add_influence(int add)
882  {
883  _influence += add;
884  return _influence;
885  }
886  // Make this point C point. Only call via amg_pointvector.
887  void make_cpoint()
888  {
889  _undecided = false;
890  _cpoint = true;
891  _influence = 0;
892  }
893  // Make this point F point. Only call via amg_pointvector.
894  void make_fpoint()
895  {
896  _undecided = false;
897  _cpoint = false;
898  _influence = 0;
899  }
900  // Switch point from F to C point. Only call via amg_pointvector.
901  void switch_ftoc() { _cpoint = true; }
902 
903  // Iterator handling for influencing and influenced points.
904  iterator begin_influencing() { return influencing_points.begin(); }
905  iterator end_influencing() { return influencing_points.end(); }
906  const_iterator begin_influencing() const { return influencing_points.begin(); }
907  const_iterator end_influencing() const { return influencing_points.end(); }
908  iterator begin_influenced() { return influenced_points.begin(); }
909  iterator end_influenced() { return influenced_points.end(); }
910  const_iterator begin_influenced() const { return influenced_points.begin(); }
911  const_iterator end_influenced() const { return influenced_points.end(); }
912  };
913 
916  struct classcomp
917  {
918  // Function returns true if l comes before r in the ordering.
919  bool operator() (amg_point* l, amg_point* r) const
920  {
921  // Map is sorted by influence number starting with the highest
922  // If influence number is the same then lowest point index comes first
923  return (l->get_influence() < r->get_influence() || (l->get_influence() == r->get_influence() && l->get_index() > r->get_index()));
924  }
925  };
926 
933  {
934  private:
935  // Type for the sorted list
936  typedef std::set<amg_point*,classcomp> ListType;
937  // Type for the vector of pointers
938  typedef std::vector<amg_point*> VectorType;
939  VectorType pointvector;
940  ListType pointlist;
941  unsigned int _size;
942  unsigned int c_points, f_points;
943 
944  public:
945  typedef VectorType::iterator iterator;
946  typedef VectorType::const_iterator const_iterator;
947 
951  amg_pointvector(unsigned int size = 0): _size(size)
952  {
953  pointvector = VectorType(size);
954  c_points = f_points = 0;
955  }
956 
957  // Construct all the points dynamically and save pointers into vector.
958  void init_points()
959  {
960  for (unsigned int i=0; i<size(); ++i)
961  pointvector[i] = new amg_point(i,size());
962  }
963  // Delete all the points.
965  {
966  for (unsigned int i=0; i<size(); ++i)
967  delete pointvector[i];
968  }
969  // 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!
970  void add_point(amg_point *point)
971  {
972  pointvector[point->get_index()] = point;
973  if (point->is_cpoint()) c_points++;
974  else if (point->is_fpoint()) f_points++;
975  }
976 
977  // Update C and F count for point *point.
978  // Necessary if C and F points were constructed outside this data structure (e.g. by parallel coarsening RS0 or RS3).
979  void update_cf(amg_point *point)
980  {
981  if (point->is_cpoint()) c_points++;
982  else if (point->is_fpoint()) f_points++;
983  }
984  // Clear the C and F point count.
985  void clear_cf() { c_points = f_points = 0; }
986 
987  // Clear both point lists.
989  {
990  for (iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
991  {
992  (*iter)->clear_influencing();
993  (*iter)->clear_influenced();
994  }
995  }
996 
997  amg_point* operator [] (unsigned int i) const { return pointvector[i]; }
998  iterator begin() { return pointvector.begin(); }
999  iterator end() { return pointvector.end(); }
1000  const_iterator begin() const { return pointvector.begin(); }
1001  const_iterator end() const { return pointvector.end(); }
1002 
1003  void resize(unsigned int size)
1004  {
1005  _size = size;
1006  pointvector = VectorType(size);
1007  }
1008  unsigned int size() const { return _size; }
1009 
1010  // Returns number of C points
1011  unsigned int get_cpoints() const { return c_points; }
1012  // Returns number of F points
1013  unsigned int get_fpoints() const { return f_points; }
1014 
1015  // Does the initial sorting of points into the list. Sorting is automatically done by the std::set data type.
1016  void sort()
1017  {
1018  for (iterator iter = begin(); iter != end(); ++iter)
1019  pointlist.insert(*iter);
1020  }
1021  // Returns the point with the highest influence measure
1023  {
1024  // No points remaining? Return NULL such that coarsening will stop.
1025  if (pointlist.size() == 0)
1026  return NULL;
1027  // If point with highest influence measure (end of the list) has measure of zero, then no further C points can be constructed. Return NULL.
1028  if ((*(--pointlist.end()))->get_influence() == 0)
1029  return NULL;
1030  // Otherwise, return the point with highest influence measure located at the end of the list.
1031  else
1032  return *(--pointlist.end());
1033  }
1034  // Add "add" to influence measure for point *point in the sorted list.
1035  void add_influence(amg_point* point, unsigned int add)
1036  {
1037  ListType::iterator iter = pointlist.find(point);
1038  // If point is not in the list then stop.
1039  if (iter == pointlist.end()) return;
1040 
1041  // Save iterator and decrement
1042  ListType::iterator iter2 = iter;
1043  iter2--;
1044 
1045  // Point has to be erased first as changing the value does not re-order the std::set
1046  pointlist.erase(iter);
1047  point->add_influence(add);
1048 
1049  // 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.
1050  pointlist.insert(iter2,point);
1051  }
1052  // Make *point to C point and remove from sorted list
1053  void make_cpoint(amg_point* point)
1054  {
1055  pointlist.erase(point);
1056  point->make_cpoint();
1057  c_points++;
1058  }
1059  // Make *point to F point and remove from sorted list
1060  void make_fpoint(amg_point* point)
1061  {
1062  pointlist.erase(point);
1063  point->make_fpoint();
1064  f_points++;
1065  }
1066  // Swich *point from F to C point
1067  void switch_ftoc(amg_point* point)
1068  {
1069  point->switch_ftoc();
1070  c_points++;
1071  f_points--;
1072  }
1073 
1074  // Build vector of indices for C point on the coarse level.
1076  {
1077  unsigned int count = 0;
1078  // Use simple counter for index creation.
1079  for (iterator iter = pointvector.begin(); iter != pointvector.end(); ++iter)
1080  {
1081  // Set index on coarse level using counter variable
1082  if ((*iter)->is_cpoint())
1083  {
1084  (*iter)->set_coarse_index(count);
1085  count++;
1086  }
1087  }
1088  }
1089 
1090  // Return information for debugging purposes
1091  template <typename MatrixType>
1092  void get_influence_matrix(MatrixType & mat) const
1093  {
1094  mat = MatrixType(size(),size());
1095  mat.clear();
1096 
1097  for (const_iterator row_iter = begin(); row_iter != end(); ++row_iter)
1098  for (amg_point::iterator col_iter = (*row_iter)->begin_influencing(); col_iter != (*row_iter)->end_influencing(); ++col_iter)
1099  mat((*row_iter)->get_index(),(*col_iter)->get_index()) = true;
1100  }
1101  template <typename VectorType>
1102  void get_influence(VectorType & vec) const
1103  {
1104  vec = VectorType(_size);
1105  vec.clear();
1106 
1107  for (const_iterator iter = begin(); iter != end(); ++iter)
1108  vec[(*iter)->get_index()] = (*iter)->get_influence();
1109  }
1110  template <typename VectorType>
1111  void get_sorting(VectorType & vec) const
1112  {
1113  vec = VectorType(pointlist.size());
1114  vec.clear();
1115  unsigned int i=0;
1116 
1117  for (ListType::const_iterator iter = pointlist.begin(); iter != pointlist.end(); ++iter)
1118  {
1119  vec[i] = (*iter)->get_index();
1120  i++;
1121  }
1122  }
1123  template <typename VectorType>
1124  void get_C(VectorType & vec) const
1125  {
1126  vec = VectorType(_size);
1127  vec.clear();
1128 
1129  for (const_iterator iter = begin(); iter != end(); ++iter)
1130  {
1131  if ((*iter)->is_cpoint())
1132  vec[(*iter)->get_index()] = true;
1133  }
1134  }
1135  template <typename VectorType>
1136  void get_F(VectorType & vec) const
1137  {
1138  vec = VectorType(_size);
1139  vec.clear();
1140 
1141  for (const_iterator iter = begin(); iter != end(); ++iter)
1142  {
1143  if ((*iter)->is_fpoint())
1144  vec[(*iter)->get_index()] = true;
1145  }
1146  }
1147  template <typename MatrixType>
1148  void get_Aggregates(MatrixType & mat) const
1149  {
1150  mat = MatrixType(_size,_size);
1151  mat.clear();
1152 
1153  for (const_iterator iter = begin(); iter != end(); ++iter)
1154  {
1155  if (!(*iter)->is_undecided())
1156  mat((*iter)->get_aggregate(),(*iter)->get_index()) = true;
1157  }
1158  }
1159  };
1160 
1164  template <typename InternalType1, typename InternalType2>
1166  {
1167  typedef typename InternalType1::value_type SparseMatrixType;
1168  typedef typename InternalType2::value_type PointVectorType;
1169 
1170  public:
1171  // Data structures on a per-processor basis.
1172  boost::numeric::ublas::vector<InternalType1> A_slice;
1173  boost::numeric::ublas::vector<InternalType2> Pointvector_slice;
1174  // Holds the offsets showing the indices for which a new slice begins.
1175  boost::numeric::ublas::vector<boost::numeric::ublas::vector<unsigned int> > Offset;
1176 
1177  unsigned int _threads;
1178  unsigned int _levels;
1179 
1180  void init(unsigned int levels, unsigned int threads = 0)
1181  {
1182  // Either use the number of threads chosen by the user or the maximum number of threads available on the processor.
1183  if (threads == 0)
1184  #ifdef _OPENMP
1185  _threads = omp_get_num_procs();
1186  #else
1187  _threads = 1;
1188  #endif
1189  else
1190  _threads = threads;
1191 
1192  _levels = levels;
1193 
1194  A_slice.resize(_threads);
1195  Pointvector_slice.resize(_threads);
1196  // Offset has _threads+1 entries to also hold the total size
1197  Offset.resize(_threads+1);
1198 
1199  for (unsigned int i=0; i<_threads; ++i)
1200  {
1201  A_slice[i].resize(_levels);
1202  Pointvector_slice[i].resize(_levels);
1203  // Offset needs one more level for the build-up of the next offset
1204  Offset[i].resize(_levels+1);
1205  }
1206  Offset[_threads].resize(_levels+1);
1207  } //init()
1208 
1209  // Slice matrix A into as many parts as threads are used.
1210  void slice (unsigned int level, InternalType1 const & A, InternalType2 const & Pointvector)
1211  {
1212  // On the finest level, build a new slicing first.
1213  if (level == 0)
1214  slice_new (level, A);
1215 
1216  // On coarser levels use the same slicing as on the finest level (Points stay together on the same thread on all levels).
1217  // This is necessary as due to interpolation and galerkin product there only exist connections between points on the same thread on coarser levels.
1218  // Note: Offset is determined in amg_coarse_rs0() after fine level was built.
1219  slice_build (level, A, Pointvector);
1220  }
1221 
1222  // Join point data structure into Pointvector
1223  void join (unsigned int level, InternalType2 & Pointvector) const
1224  {
1225  typedef typename InternalType2::value_type PointVectorType;
1226 
1227  // Reset index offset of all points and update overall C and F point count
1228  Pointvector[level].clear_cf();
1229  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
1230  {
1231  (*iter)->set_offset(0);
1232  Pointvector[level].update_cf(*iter);
1233  }
1234  }
1235 
1236  private:
1241  void slice_new (unsigned int level, InternalType1 const & A)
1242  {
1243  typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
1244  typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
1245 
1246  // Determine index offset of all the slices (index of A[level] when the respective slice starts).
1247  #ifdef _OPENMP
1248  #pragma omp parallel for
1249  #endif
1250  for (unsigned int i=0; i<=_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] = A[level].size1();
1255  else Offset[i][level] = 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 _OPENMP
1273  #pragma omp parallel for private (x,y,point)
1274  #endif
1275  for (unsigned int i=0; i<_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 = 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 = 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 = col_iter.index2();
1304  }
1305 
1306  ++row_iter;
1307  x = 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  unsigned int x,y,z;
1326  ScalarType prod;
1327  RES = SparseMatrixType(A.size1(), B.size2());
1328  RES.clear();
1329 
1330  #ifdef _OPENMP
1331  #pragma omp parallel for private (x,y,z,prod) shared (A,B,RES)
1332  #endif
1333  for (x=0; x<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 = 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 = 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  unsigned int x,y1,y2,z;
1367  RES = SparseMatrixType(P.size2(), P.size2());
1368  RES.clear();
1369 
1370  #ifdef _OPENMP
1371  #pragma omp parallel for private (x,y1,y2,z,row) shared (A,P,RES)
1372  #endif
1373  for (x=0; x<P.size2(); ++x)
1374  {
1375  row = amg_sparsevector<ScalarType>(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 = 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 = 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 = col_iter3.index2();
1400  RES.add (x, z, *col_iter3 * *iter);
1401  }
1402  }
1403  }
1404 
1405  #ifdef 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 (abs((ScalarType)RAP(x,y) - (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