ViennaCL - The Vienna Computing Library  1.2.0
amg_coarse.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
2 #define VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
3 
4 /* =========================================================================
5  Copyright (c) 2010-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
24 #include <cmath>
25 #include "viennacl/linalg/amg.hpp"
26 
27 #include <map>
28 #ifdef _OPENMP
29 #include <omp.h>
30 #endif
31 
32 #include "amg_debug.hpp"
33 
34 namespace viennacl
35 {
36  namespace linalg
37  {
38  namespace detail
39  {
40  namespace amg
41  {
42 
50  template <typename InternalType1, typename InternalType2, typename InternalType3>
51  void amg_coarse(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing, amg_tag & tag)
52  {
53  switch (tag.get_coarse())
54  {
55  case VIENNACL_AMG_COARSE_RS: amg_coarse_classic (level, A, Pointvector, tag); break;
56  case VIENNACL_AMG_COARSE_ONEPASS: amg_coarse_classic_onepass (level, A, Pointvector, tag); break;
57  case VIENNACL_AMG_COARSE_RS0: amg_coarse_rs0 (level, A, Pointvector, Slicing, tag); break;
58  case VIENNACL_AMG_COARSE_RS3: amg_coarse_rs3 (level, A, Pointvector, Slicing, tag); break;
59  case VIENNACL_AMG_COARSE_AG: amg_coarse_ag (level, A, Pointvector, tag); break;
60  }
61  }
62 
69  template <typename InternalType1, typename InternalType2>
70  void amg_influence(unsigned int level, InternalType1 const & A, InternalType2 & Pointvector, amg_tag & tag)
71  {
72  typedef typename InternalType1::value_type SparseMatrixType;
73  typedef typename InternalType2::value_type PointVectorType;
74  typedef typename SparseMatrixType::value_type ScalarType;
75  typedef typename SparseMatrixType::value_type ScalarType;
76  typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
77  typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
78 
79  ScalarType max;
80  int diag_sign;
81  //unsigned int i;
82 
83 #ifdef _OPENMP
84  #pragma omp parallel for private (max,diag_sign) shared (A,Pointvector)
85 #endif
86  for (unsigned int i=0; i<A[level].size1(); ++i)
87  {
88  diag_sign = 1;
89  if (A[level](i,i) < 0)
90  diag_sign = -1;
91 
92  ConstRowIterator row_iter = A[level].begin1();
93  row_iter += i;
94  // Find greatest non-diagonal negative value (positive if diagonal is negative) in row
95  max = 0;
96  for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
97  {
98  if (i == (unsigned int) col_iter.index2()) continue;
99  if (diag_sign == 1)
100  if (max > *col_iter) max = *col_iter;
101  if (diag_sign == -1)
102  if (max < *col_iter) max = *col_iter;
103  }
104 
105  // If maximum is 0 then the row is independent of the others
106  if (max == 0)
107  continue;
108 
109  // Find all points that strongly influence current point (Yang, p.5)
110  for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
111  {
112  unsigned int j = col_iter.index2();
113  if (i == j) continue;
114  if (diag_sign * (-*col_iter) >= tag.get_threshold() * (diag_sign * (-max)))
115  {
116  // Strong influence from j to i found, save information
117  Pointvector[level][i]->add_influencing_point(Pointvector[level][j]);
118  }
119  }
120  }
121 
122  #ifdef DEBUG
123  std::cout << "Influence Matrix: " << std::endl;
124  boost::numeric::ublas::matrix<bool> mat;
125  Pointvector[level].get_influence_matrix(mat);
126  printmatrix (mat);
127  #endif
128 
129  // Save influenced points
130  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
131  {
132  for (typename amg_point::iterator iter2 = (*iter)->begin_influencing(); iter2 != (*iter)->end_influencing(); ++iter2)
133  {
134  (*iter2)->add_influenced_point(*iter);
135  }
136  }
137 
138  #ifdef DEBUG
139  std::cout << "Influence Measures: " << std::endl;
140  boost::numeric::ublas::vector<unsigned int> temp;
141  Pointvector[level].get_influence(temp);
142  printvector (temp);
143  std::cout << "Point Sorting: " << std::endl;
144  Pointvector[level].get_sorting(temp);
145  printvector (temp);
146  #endif
147  }
148 
155  template <typename InternalType1, typename InternalType2>
156  void amg_coarse_classic_onepass(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, amg_tag & tag)
157  {
158  typedef typename InternalType1::value_type SparseMatrixType;
159  typedef typename InternalType2::value_type PointVectorType;
160  typedef typename SparseMatrixType::value_type ScalarType;
161 
162  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
163  typedef typename SparseMatrixType::iterator2 InternalColIterator;
164 
165  amg_point* c_point, *point1, *point2;
166  unsigned int i;
167 
168  // Check and save all strong influences
169  amg_influence (level, A, Pointvector, tag);
170 
171  // Traverse through points and calculate initial influence measure
172 #ifdef _OPENMP
173  #pragma omp parallel for private (i) shared (Pointvector)
174 #endif
175  for (i=0; i<Pointvector[level].size(); ++i)
176  Pointvector[level][i]->calc_influence();
177 
178  // Do initial sorting
179  Pointvector[level].sort();
180 
181  // Get undecided point with highest influence measure
182  while ((c_point = Pointvector[level].get_nextpoint()) != NULL)
183  {
184  // Make this point C point
185  Pointvector[level].make_cpoint(c_point);
186 
187  // All strongly influenced points become F points
188  for (typename amg_point::iterator iter = c_point->begin_influenced(); iter != c_point->end_influenced(); ++iter)
189  {
190  point1 = *iter;
191  // Found strong influence from C point (c_point influences point1), check whether point is still undecided, otherwise skip
192  if (!point1->is_undecided()) continue;
193  // Make this point F point if it is still undecided point
194  Pointvector[level].make_fpoint(point1);
195 
196  // Add +1 to influence measure for all undecided points that strongly influence new F point
197  for (typename amg_point::iterator iter2 = point1->begin_influencing(); iter2 != point1->end_influencing(); ++iter2)
198  {
199  point2 = *iter2;
200  // Found strong influence to F point (point2 influences point1)
201  if (point2->is_undecided())
202  Pointvector[level].add_influence(point2,1);
203  }
204  }
205  }
206 
207  // If a point is neither C nor F point but is nevertheless influenced by other points make it F point
208  // (this situation can happen when this point does not influence other points and the points that influence this point became F points already)
209  /*#pragma omp parallel for private (i,point1)
210  for (i=0; i<Pointvector[level].size(); ++i)
211  {
212  point1 = Pointvector[level][i];
213  if (point1->is_undecided())
214  {
215  // Undecided point found. Check whether it is influenced by other point and if so: Make it F point.
216  if (point1->number_influencing() > 0)
217  {
218  #pragma omp critical
219  Pointvector[level].make_fpoint(point1);
220  }
221  }
222  }*/
223 
224  #if defined (DEBUG)// or defined (DEBUGBENCH)
225  unsigned int c_points = Pointvector[level].get_cpoints();
226  unsigned int f_points = Pointvector[level].get_fpoints();
227  std::cout << "1st pass: Level " << level << ": ";
228  std::cout << "No of C points = " << c_points << ", ";
229  std::cout << "No of F points = " << f_points << std::endl;
230  #endif
231 
232  #ifdef DEBUG
233  std::cout << "Coarse Points:" << std::endl;
234  boost::numeric::ublas::vector<bool> C;
235  Pointvector[level].get_C(C);
236  printvector (C);
237  std::cout << "Fine Points:" << std::endl;
238  boost::numeric::ublas::vector<bool> F;
239  Pointvector[level].get_F(F);
240  printvector (F);
241  #endif
242  }
243 
250  template <typename InternalType1, typename InternalType2>
251  void amg_coarse_classic(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, amg_tag & tag)
252  {
253  typedef typename InternalType1::value_type SparseMatrixType;
254  typedef typename InternalType2::value_type PointVectorType;
255  typedef typename SparseMatrixType::value_type ScalarType;
256 
257  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
258  typedef typename SparseMatrixType::iterator2 InternalColIterator;
259 
260  bool add_C;
261  amg_point *c_point, *point1, *point2;
262 
263  // Use one-pass-coarsening as first pass.
264  amg_coarse_classic_onepass(level, A, Pointvector, tag);
265 
266  // 2nd pass: Add more C points if F-F connection does not have a common C point.
267  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
268  {
269  point1 = *iter;
270  // If point is F point, check for strong connections.
271  if (point1->is_fpoint())
272  {
273  // Check for strong connections from influencing and influenced points.
274  amg_point::iterator iter2 = point1->begin_influencing();
275  amg_point::iterator iter3 = point1->begin_influenced();
276 
277  // Iterate over both lists at once. This makes sure that points are no checked twice when influence relation is symmetric (which is often the case).
278  // Note: Only works because influencing and influenced lists are sorted by point-index.
279  while(iter2 != point1->end_influencing() || iter3 != point1->end_influenced())
280  {
281  if (iter2 == point1->end_influencing())
282  {
283  point2 = *iter3;
284  ++iter3;
285  }
286  else if (iter3 == point1->end_influenced())
287  {
288  point2 = *iter2;
289  ++iter2;
290  }
291  else
292  {
293  if ((*iter2)->get_index() == (*iter3)->get_index())
294  {
295  point2 = *iter2;
296  ++iter2;
297  ++iter3;
298  }
299  else if ((*iter2)->get_index() < (*iter3)->get_index())
300  {
301  point2 = *iter2;
302  ++iter2;
303  }
304  else
305  {
306  point2 = *iter3;
307  ++iter3;
308  }
309  }
310  // Only check points with higher index as points with lower index have been checked already.
311  if (point2->get_index() < point1->get_index())
312  continue;
313 
314  // If there is a strong connection then it has to either be a C point or a F point with common C point.
315  // C point? Then skip as everything is ok.
316  if (point2->is_cpoint())
317  continue;
318  // F point? Then check whether F points point1 and point2 have a common C point.
319  if (point2->is_fpoint())
320  {
321  add_C = true;
322  // C point is common for two F points if they are both strongly influenced by that C point.
323  // Compare strong influences for point1 and point2.
324  for (amg_point::iterator iter3 = point1->begin_influencing(); iter3 != point1 -> end_influencing(); ++iter3)
325  {
326  c_point = *iter3;
327  // Stop search when strong common influence is found via c_point.
328  if (c_point->is_cpoint())
329  {
330  if (point2->is_influencing(c_point))
331  {
332  add_C = false;
333  break;
334  }
335  }
336  }
337  // No common C point found? Then make second F point to C point.
338  if (add_C == true)
339  Pointvector[level].switch_ftoc(point2);
340  }
341  }
342  }
343  }
344 
345  #ifdef DEBUG
346  std::cout << "After 2nd pass:" << std::endl;
347  std::cout << "Coarse Points:" << std::endl;
348  boost::numeric::ublas::vector<bool> C;
349  Pointvector[level].get_C(C);
350  printvector (C);
351  std::cout << "Fine Points:" << std::endl;
352  boost::numeric::ublas::vector<bool> F;
353  Pointvector[level].get_F(F);
354  printvector (F);
355  #endif
356 
357  #ifdef DEBUG
358 #ifdef _OPENMP
359  #pragma omp critical
360 #endif
361  {
362  std::cout << "No C and no F point: ";
363  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
364  if ((*iter)->is_undecided())
365  std::cout << (*iter)->get_index() << " ";
366  std::cout << std::endl;
367  }
368  #endif
369  }
370 
378  template <typename InternalType1, typename InternalType2, typename InternalType3>
379  void amg_coarse_rs0(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing, amg_tag & tag)
380  {
381  typedef typename InternalType1::value_type SparseMatrixType;
382  typedef typename InternalType2::value_type PointVectorType;
383  typedef typename SparseMatrixType::value_type ScalarType;
384 
385  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
386  typedef typename SparseMatrixType::iterator2 InternalColIterator;
387 
388  unsigned int total_points;
389 
390  // Slice matrix into parts such that points are distributed among threads
391  Slicing.slice(level, A, Pointvector);
392 
393  // Run classical coarsening in parallel
394  total_points = 0;
395 #ifdef _OPENMP
396  #pragma omp parallel for shared (total_points,Slicing,level)
397 #endif
398  for (unsigned int i=0; i<Slicing._threads; ++i)
399  {
400  amg_coarse_classic(level,Slicing.A_slice[i],Slicing.Pointvector_slice[i],tag);
401 
402  // Save C points (using Slicing.Offset on the next level as temporary memory)
403  // Note: Number of C points for point i is saved in i+1!! (makes it easier later to compute offset)
404  Slicing.Offset[i+1][level+1] = Slicing.Pointvector_slice[i][level].get_cpoints();
405 #ifdef _OPENMP
406  #pragma omp critical
407 #endif
408  total_points += Slicing.Pointvector_slice[i][level].get_cpoints();
409  }
410 
411  // If no coarser level can be found on any level then resume and coarsening will stop in amg_coarse()
412  if (total_points != 0)
413  {
414 #ifdef _OPENMP
415  #pragma omp parallel for shared (Slicing)
416 #endif
417  for (unsigned int i=0; i<Slicing._threads; ++i)
418  {
419  // If no higher coarse level can be found on slice i (saved in Slicing.Offset[i+1][level+1]) then pull C point(s) to the next level
420  if (Slicing.Offset[i+1][level+1] == 0)
421  {
422  // All points become C points
423  for (unsigned int j=0; j<Slicing.A_slice[i][level].size1(); ++j)
424  Slicing.Pointvector_slice[i][level].make_cpoint(Slicing.Pointvector_slice[i][level][j]);
425  Slicing.Offset[i+1][level+1] = Slicing.A_slice[i][level].size1();
426  }
427  }
428 
429  // Build slicing offset from number of C points (offset = total sum of C points on threads with lower number)
430  for (unsigned int i=2; i<=Slicing._threads; ++i)
431  Slicing.Offset[i][level+1] += Slicing.Offset[i-1][level+1];
432 
433  // Join C and F points
434  Slicing.join(level, Pointvector);
435  }
436 
437  // Calculate global influence measures for interpolation and/or RS3.
438  amg_influence(level, A, Pointvector, tag);
439 
440  #if defined(DEBUG)// or defined (DEBUGBENCH)
441  for (unsigned int i=0; i<Slicing._threads; ++i)
442  {
443  unsigned int c_points = Slicing.Pointvector_slice[i][level].get_cpoints();
444  unsigned int f_points = Slicing.Pointvector_slice[i][level].get_fpoints();
445  std::cout << "Thread " << i << ": ";
446  std::cout << "No of C points = " << c_points << ", ";
447  std::cout << "No of F points = " << f_points << std::endl;
448  }
449  #endif
450  }
451 
459  template <typename InternalType1, typename InternalType2, typename InternalType3>
460  void amg_coarse_rs3(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing, amg_tag & tag)
461  {
462  typedef typename InternalType1::value_type SparseMatrixType;
463  typedef typename InternalType2::value_type PointVectorType;
464  typedef typename SparseMatrixType::value_type ScalarType;
465 
466  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
467  typedef typename SparseMatrixType::iterator2 InternalColIterator;
468 
469  amg_point *c_point, *point1, *point2;
470  bool add_C;
471  unsigned int i, j;
472 
473  // Run RS0 first (parallel).
474  amg_coarse_rs0(level, A, Pointvector, Slicing, tag);
475 
476  // Save slicing offset
477  boost::numeric::ublas::vector<unsigned int> Offset = boost::numeric::ublas::vector<unsigned int> (Slicing.Offset.size());
478  for (i=0; i<Slicing.Offset.size(); ++i)
479  Offset[i] = Slicing.Offset[i][level];
480 
481  // Correct the coarsening with a third pass: Don't allow strong F-F connections without common C point
482  for (i=0; i<Slicing._threads; ++i)
483  {
484  //for (j=Slicing.Offset[i][level]; j<Slicing.Offset[i+1][level]; ++j)
485  for (j=Offset[i]; j<Offset[i+1]; ++j)
486  {
487  point1 = Pointvector[level][j];
488  // If point is F point, check for strong connections.
489  if (point1->is_fpoint())
490  {
491  // Check for strong connections from influencing and influenced points.
492  amg_point::iterator iter2 = point1->begin_influencing();
493  amg_point::iterator iter3 = point1->begin_influenced();
494 
495  // Iterate over both lists at once. This makes sure that points are no checked twice when influence relation is symmetric (which is often the case).
496  // Note: Only works because influencing and influenced lists are sorted by point-index.
497  while(iter2 != point1->end_influencing() || iter3 != point1->end_influenced())
498  {
499  if (iter2 == point1->end_influencing())
500  {
501  point2 = *iter3;
502  ++iter3;
503  }
504  else if (iter3 == point1->end_influenced())
505  {
506  point2 = *iter2;
507  ++iter2;
508  }
509  else
510  {
511  if ((*iter2)->get_index() == (*iter3)->get_index())
512  {
513  point2 = *iter2;
514  ++iter2;
515  ++iter3;
516  }
517  else if ((*iter2)->get_index() < (*iter3)->get_index())
518  {
519  point2 = *iter2;
520  ++iter2;
521  }
522  else
523  {
524  point2 = *iter3;
525  ++iter3;
526  }
527  }
528 
529  // Only check points with higher index as points with lower index have been checked already.
530  if (point2->get_index() < point1->get_index())
531  continue;
532 
533  // Only check points that are outside the slicing boundaries (interior F-F connections have already been checked in second pass)
534  //if (point2->get_index() >= Slicing.Offset[i][level] || point2->get_index() < Slicing.Offset[i+1][level])
535  if (point2->get_index() >= Offset[i] && point2->get_index() < Offset[i+1])
536  continue;
537 
538  // If there is a strong connection then it has to either be a C point or a F point with common C point.
539  // C point? Then skip as everything is ok.
540  if (point2->is_cpoint())
541  continue;
542  // F point? Then check whether F points point1 and point2 have a common C point.
543  if (point2->is_fpoint())
544  {
545  add_C = true;
546  // C point is common for two F points if they are both strongly influenced by that C point.
547  // Compare strong influences for point1 and point2.
548  for (amg_point::iterator iter3 = point1->begin_influencing(); iter3 != point1 -> end_influencing(); ++iter3)
549  {
550  c_point = *iter3;
551  // Stop search when strong common influence is found via c_point.
552  if (c_point->is_cpoint())
553  {
554  if (point2->is_influencing(c_point))
555  {
556  add_C = false;
557  break;
558  }
559  }
560  }
561  // No common C point found? Then make second F point to C point.
562  if (add_C == true)
563  {
564  Pointvector[level].switch_ftoc(point2);
565  // Add +1 to offsets as one C point has been added.
566  for (unsigned int j=i+1; j<=Slicing._threads; ++j)
567  Slicing.Offset[j][level+1]++;
568  }
569  }
570  }
571  }
572  }
573  }
574 
575  #ifdef DEBUG
576  std::cout << "After 3rd pass:" << std::endl;
577  std::cout << "Coarse Points:" << std::endl;
578  boost::numeric::ublas::vector<bool> C;
579  Pointvector[level].get_C(C);
580  printvector (C);
581  std::cout << "Fine Points:" << std::endl;
582  boost::numeric::ublas::vector<bool> F;
583  Pointvector[level].get_F(F);
584  printvector (F);
585  #endif
586 
587  #ifdef DEBUG
588  unsigned int i;
589 #ifdef _OPENMP
590  #pragma omp critical
591 #endif
592  {
593  std::cout << "No C and no F point: ";
594  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
595  if ((*iter)->is_undecided())
596  std::cout << i << " ";
597  std::cout << std::endl;
598  }
599  #endif
600  }
601 
609  template <typename InternalType1, typename InternalType2>
610  void amg_coarse_ag(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, amg_tag & tag)
611  {
612  typedef typename InternalType1::value_type SparseMatrixType;
613  typedef typename InternalType2::value_type PointVectorType;
614  typedef typename SparseMatrixType::value_type ScalarType;
615 
616  typedef typename SparseMatrixType::iterator1 InternalRowIterator;
617  typedef typename SparseMatrixType::iterator2 InternalColIterator;
618 
619  unsigned int x,y;
620  ScalarType diag;
621  amg_point *pointx, *pointy;
622 
623  // Cannot determine aggregates if size == 1 as then a new aggregate would always consist of this point (infinite loop)
624  if (A[level].size1() == 1) return;
625 
626  // SA algorithm (Vanek et al. p.6)
627  // Build neighborhoods
628 #ifdef _OPENMP
629  #pragma omp parallel for private (x,y,diag) shared (A)
630 #endif
631  for (x=0; x<A[level].size1(); ++x)
632  {
633  InternalRowIterator row_iter = A[level].begin1();
634  row_iter += x;
635  diag = A[level](x,x);
636  for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
637  {
638  y = col_iter.index2();
639  if (y == x || (std::abs(*col_iter) >= tag.get_threshold()*pow(0.5,level-1) * sqrt(std::abs(diag*A[level](y,y)))))
640  {
641  // Neighborhood x includes point y
642  Pointvector[level][x]->add_influencing_point(Pointvector[level][y]);
643  }
644  }
645  }
646 
647  #ifdef DEBUG
648  std::cout << "Neighborhoods:" << std::endl;
649  boost::numeric::ublas::matrix<bool> mat;
650  Pointvector[level].get_influence_matrix(mat);
651  printmatrix (mat);
652  #endif
653 
654  // Build aggregates from neighborhoods
655  for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
656  {
657  pointx = (*iter);
658 
659  if (pointx->is_undecided())
660  {
661  // Make center of aggregate to C point and include it to aggregate x.
662  Pointvector[level].make_cpoint(pointx);
663  pointx->set_aggregate (pointx->get_index());
664  for (amg_point::iterator iter2 = pointx->begin_influencing(); iter2 != pointx->end_influencing(); ++iter2)
665  {
666  pointy = (*iter2);
667 
668  if (pointy->is_undecided())
669  {
670  // Make neighbor y to F point and include it to aggregate x.
671  Pointvector[level].make_fpoint(pointy);
672  pointy->set_aggregate (pointx->get_index());
673  }
674  }
675  }
676  }
677 
678  #ifdef DEBUG
679  std::cout << "After aggregation:" << std::endl;
680  std::cout << "Coarse Points:" << std::endl;
681  boost::numeric::ublas::vector<bool> C;
682  Pointvector[level].get_C(C);
683  printvector (C);
684  std::cout << "Fine Points:" << std::endl;
685  boost::numeric::ublas::vector<bool> F;
686  Pointvector[level].get_F(F);
687  printvector (F);
688  std::cout << "Aggregates:" << std::endl;
689  printvector (Aggregates[level]);
690  #endif
691  }
692  } //namespace amg
693  }
694  }
695 }
696 
697 #endif