Reference documentation for deal.II version 8.4.1
mg_smoother.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2016 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__mg_smoother_h
17 #define dealii__mg_smoother_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/smartpointer.h>
22 #include <deal.II/lac/pointer_matrix.h>
23 #include <deal.II/lac/vector_memory.h>
24 #include <deal.II/multigrid/mg_base.h>
25 #include <deal.II/base/mg_level_object.h>
26 #include <vector>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 /*
31  * MGSmootherBase is defined in mg_base.h
32  */
33 
36 
45 template <typename VectorType>
46 class MGSmoother : public MGSmootherBase<VectorType>
47 {
48 public:
52  MGSmoother(const unsigned int steps = 1,
53  const bool variable = false,
54  const bool symmetric = false,
55  const bool transpose = false);
56 
60  void set_steps (const unsigned int);
61 
65  void set_variable (const bool);
66 
70  void set_symmetric (const bool);
71 
76  void set_transpose (const bool);
77 
82  void set_debug (const unsigned int level);
83 
84 protected:
92 
97  unsigned int steps;
98 
103  bool variable;
104 
109  bool symmetric;
110 
115  bool transpose;
116 
120  unsigned int debug;
121 };
122 
123 
132 template <typename VectorType>
133 class MGSmootherIdentity : public MGSmootherBase<VectorType>
134 {
135 public:
141  virtual void smooth (const unsigned int level,
142  VectorType &u,
143  const VectorType &rhs) const;
144  virtual void clear ();
145 };
146 
147 
148 namespace mg
149 {
179  template<class RelaxationType, typename VectorType>
180  class SmootherRelaxation : public MGLevelObject<RelaxationType>, public MGSmoother<VectorType>
181  {
182  public:
186  SmootherRelaxation(const unsigned int steps = 1,
187  const bool variable = false,
188  const bool symmetric = false,
189  const bool transpose = false);
190 
199  template <typename MatrixType2>
200  void initialize (const MGLevelObject<MatrixType2> &matrices,
201  const typename RelaxationType::AdditionalData &additional_data
202  = typename RelaxationType::AdditionalData());
203 
211  template <typename MatrixType2, class DATA>
212  void initialize (const MGLevelObject<MatrixType2> &matrices,
213  const MGLevelObject<DATA> &additional_data);
214 
218  void clear ();
219 
223  virtual void smooth (const unsigned int level,
224  VectorType &u,
225  const VectorType &rhs) const;
226 
230  std::size_t memory_consumption () const;
231  };
232 }
233 
265 template<typename MatrixType, class RelaxationType, typename VectorType>
266 class MGSmootherRelaxation : public MGSmoother<VectorType>
267 {
268 public:
272  MGSmootherRelaxation(const unsigned int steps = 1,
273  const bool variable = false,
274  const bool symmetric = false,
275  const bool transpose = false);
276 
285  template <typename MatrixType2>
287  const typename RelaxationType::AdditionalData &additional_data
288  = typename RelaxationType::AdditionalData());
289 
298  template <typename MatrixType2, class DATA>
299  void initialize (const MGLevelObject<MatrixType2> &matrices,
300  const MGLevelObject<DATA> &additional_data);
301 
311  template <typename MatrixType2, class DATA>
312  void initialize (const MGLevelObject<MatrixType2> &matrices,
313  const DATA &additional_data,
314  const unsigned int block_row,
315  const unsigned int block_col);
316 
326  template <typename MatrixType2, class DATA>
327  void initialize (const MGLevelObject<MatrixType2> &matrices,
328  const MGLevelObject<DATA> &additional_data,
329  const unsigned int block_row,
330  const unsigned int block_col);
331 
335  void clear ();
336 
340  virtual void smooth (const unsigned int level,
341  VectorType &u,
342  const VectorType &rhs) const;
343 
348 
352  std::size_t memory_consumption () const;
353 
354 
355 private:
360 
361 };
362 
363 
364 
395 template<typename MatrixType, typename PreconditionerType, typename VectorType>
396 class MGSmootherPrecondition : public MGSmoother<VectorType>
397 {
398 public:
402  MGSmootherPrecondition(const unsigned int steps = 1,
403  const bool variable = false,
404  const bool symmetric = false,
405  const bool transpose = false);
406 
416  template <typename MatrixType2>
418  const typename PreconditionerType::AdditionalData &additional_data = typename PreconditionerType::AdditionalData());
419 
429  template <typename MatrixType2, class DATA>
430  void initialize (const MGLevelObject<MatrixType2> &matrices,
431  const MGLevelObject<DATA> &additional_data);
432 
443  template <typename MatrixType2, class DATA>
444  void initialize (const MGLevelObject<MatrixType2> &matrices,
445  const DATA &additional_data,
446  const unsigned int block_row,
447  const unsigned int block_col);
448 
459  template <typename MatrixType2, class DATA>
460  void initialize (const MGLevelObject<MatrixType2> &matrices,
461  const MGLevelObject<DATA> &additional_data,
462  const unsigned int block_row,
463  const unsigned int block_col);
464 
468  void clear ();
469 
473  virtual void smooth (const unsigned int level,
474  VectorType &u,
475  const VectorType &rhs) const;
476 
481 
485  std::size_t memory_consumption () const;
486 
487 
488 private:
493 
494 };
495 
498 /* ------------------------------- Inline functions -------------------------- */
499 
500 #ifndef DOXYGEN
501 
502 template <typename VectorType>
503 inline void
504 MGSmootherIdentity<VectorType>::smooth (const unsigned int,
505  VectorType &,
506  const VectorType &) const
507 {}
508 
509 template <typename VectorType>
510 inline void
512 {}
513 
514 //---------------------------------------------------------------------------
515 
516 template <typename VectorType>
517 inline
518 MGSmoother<VectorType>::MGSmoother (const unsigned int steps,
519  const bool variable,
520  const bool symmetric,
521  const bool transpose)
522  :
523  steps(steps),
524  variable(variable),
525  symmetric(symmetric),
526  transpose(transpose),
527  debug(0)
528 {}
529 
530 
531 template <typename VectorType>
532 inline void
533 MGSmoother<VectorType>::set_steps (const unsigned int s)
534 {
535  steps = s;
536 }
537 
538 
539 template <typename VectorType>
540 inline void
541 MGSmoother<VectorType>::set_debug (const unsigned int s)
542 {
543  debug = s;
544 }
545 
546 
547 template <typename VectorType>
548 inline void
549 MGSmoother<VectorType>::set_variable (const bool flag)
550 {
551  variable = flag;
552 }
553 
554 
555 template <typename VectorType>
556 inline void
558 {
559  symmetric = flag;
560 }
561 
562 
563 template <typename VectorType>
564 inline void
566 {
567  transpose = flag;
568 }
569 
570 //----------------------------------------------------------------------//
571 
572 namespace mg
573 {
574  template <class RelaxationType, typename VectorType>
575  inline
577  (const unsigned int steps,
578  const bool variable,
579  const bool symmetric,
580  const bool transpose)
581  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
582  {}
583 
584 
585  template <class RelaxationType, typename VectorType>
586  inline void
587  SmootherRelaxation<RelaxationType, VectorType>::clear ()
588  {
590  }
591 
592 
593  template <class RelaxationType, typename VectorType>
594  template <typename MatrixType2>
595  inline void
596  SmootherRelaxation<RelaxationType, VectorType>::initialize
597  (const MGLevelObject<MatrixType2> &m,
598  const typename RelaxationType::AdditionalData &data)
599  {
600  const unsigned int min = m.min_level();
601  const unsigned int max = m.max_level();
602 
603  this->resize(min, max);
604 
605  for (unsigned int i=min; i<=max; ++i)
606  (*this)[i].initialize(m[i], data);
607  }
608 
609 
610  template <class RelaxationType, typename VectorType>
611  template <typename MatrixType2, class DATA>
612  inline void
613  SmootherRelaxation<RelaxationType, VectorType>::initialize
614  (const MGLevelObject<MatrixType2> &m,
615  const MGLevelObject<DATA> &data)
616  {
617  const unsigned int min = std::max(m.min_level(), data.min_level());
618  const unsigned int max = std::min(m.max_level(), data.max_level());
619 
620  this->resize(min, max);
621 
622  for (unsigned int i=min; i<=max; ++i)
623  (*this)[i].initialize(m[i], data[i]);
624  }
625 
626 
627  template <class RelaxationType, typename VectorType>
628  inline void
629  SmootherRelaxation<RelaxationType, VectorType>::smooth (const unsigned int level,
630  VectorType &u,
631  const VectorType &rhs) const
632  {
633  unsigned int maxlevel = this->max_level();
634  unsigned int steps2 = this->steps;
635 
636  if (this->variable)
637  steps2 *= (1<<(maxlevel-level));
638 
639  bool T = this->transpose;
640  if (this->symmetric && (steps2 % 2 == 0))
641  T = false;
642  if (this->debug > 0)
643  deallog << 'S' << level << ' ';
644 
645  for (unsigned int i=0; i<steps2; ++i)
646  {
647  if (T)
648  (*this)[level].Tstep(u, rhs);
649  else
650  (*this)[level].step(u, rhs);
651  if (this->symmetric)
652  T = !T;
653  }
654  }
655 
656 
657  template <class RelaxationType, typename VectorType>
658  inline
659  std::size_t
660  SmootherRelaxation<RelaxationType, VectorType>::
661  memory_consumption () const
662  {
663  return sizeof(*this)
666  + this->vector_memory.memory_consumption();
667  }
668 }
669 
670 
671 //----------------------------------------------------------------------//
672 
673 template <typename MatrixType, class RelaxationType, typename VectorType>
674 inline
676 (const unsigned int steps,
677  const bool variable,
678  const bool symmetric,
679  const bool transpose)
680  :
681  MGSmoother<VectorType>(steps, variable, symmetric, transpose)
682 {}
683 
684 
685 
686 template <typename MatrixType, class RelaxationType, typename VectorType>
687 inline void
689 {
690  smoothers.clear();
691 
692  unsigned int i=matrices.min_level(),
693  max_level=matrices.max_level();
694  for (; i<=max_level; ++i)
695  matrices[i]=0;
696 }
697 
698 
699 template <typename MatrixType, class RelaxationType, typename VectorType>
700 template <typename MatrixType2>
701 inline void
704  const typename RelaxationType::AdditionalData &data)
705 {
706  const unsigned int min = m.min_level();
707  const unsigned int max = m.max_level();
708 
709  matrices.resize(min, max);
710  smoothers.resize(min, max);
711 
712  for (unsigned int i=min; i<=max; ++i)
713  {
714  matrices[i] = &m[i];
715  smoothers[i].initialize(m[i], data);
716  }
717 }
718 
719 template <typename MatrixType, class RelaxationType, typename VectorType>
720 template <typename MatrixType2, class DATA>
721 inline void
724  const MGLevelObject<DATA> &data)
725 {
726  const unsigned int min = m.min_level();
727  const unsigned int max = m.max_level();
728 
729  Assert (data.min_level() == min,
730  ExcDimensionMismatch(data.min_level(), min));
731  Assert (data.max_level() == max,
732  ExcDimensionMismatch(data.max_level(), max));
733 
734  matrices.resize(min, max);
735  smoothers.resize(min, max);
736 
737  for (unsigned int i=min; i<=max; ++i)
738  {
739  matrices[i] = &m[i];
740  smoothers[i].initialize(m[i], data[i]);
741  }
742 }
743 
744 template <typename MatrixType, class RelaxationType, typename VectorType>
745 template <typename MatrixType2, class DATA>
746 inline void
749  const DATA &data,
750  const unsigned int row,
751  const unsigned int col)
752 {
753  const unsigned int min = m.min_level();
754  const unsigned int max = m.max_level();
755 
756  matrices.resize(min, max);
757  smoothers.resize(min, max);
758 
759  for (unsigned int i=min; i<=max; ++i)
760  {
761  matrices[i] = &(m[i].block(row, col));
762  smoothers[i].initialize(m[i].block(row, col), data);
763  }
764 }
765 
766 template <typename MatrixType, class RelaxationType, typename VectorType>
767 template <typename MatrixType2, class DATA>
768 inline void
771  const MGLevelObject<DATA> &data,
772  const unsigned int row,
773  const unsigned int col)
774 {
775  const unsigned int min = m.min_level();
776  const unsigned int max = m.max_level();
777 
778  Assert (data.min_level() == min,
779  ExcDimensionMismatch(data.min_level(), min));
780  Assert (data.max_level() == max,
781  ExcDimensionMismatch(data.max_level(), max));
782 
783  matrices.resize(min, max);
784  smoothers.resize(min, max);
785 
786  for (unsigned int i=min; i<=max; ++i)
787  {
788  matrices[i] = &(m[i].block(row, col));
789  smoothers[i].initialize(m[i].block(row, col), data[i]);
790  }
791 }
792 
793 
794 template <typename MatrixType, class RelaxationType, typename VectorType>
795 inline void
797  VectorType &u,
798  const VectorType &rhs) const
799 {
800  unsigned int maxlevel = smoothers.max_level();
801  unsigned int steps2 = this->steps;
802 
803  if (this->variable)
804  steps2 *= (1<<(maxlevel-level));
805 
806  bool T = this->transpose;
807  if (this->symmetric && (steps2 % 2 == 0))
808  T = false;
809  if (this->debug > 0)
810  deallog << 'S' << level << ' ';
811 
812  for (unsigned int i=0; i<steps2; ++i)
813  {
814  if (T)
815  smoothers[level].Tstep(u, rhs);
816  else
817  smoothers[level].step(u, rhs);
818  if (this->symmetric)
819  T = !T;
820  }
821 }
822 
823 
824 
825 template <typename MatrixType, class RelaxationType, typename VectorType>
826 inline
827 std::size_t
829 memory_consumption () const
830 {
831  return sizeof(*this)
832  + matrices.memory_consumption()
833  + smoothers.memory_consumption()
834  + this->vector_memory.memory_consumption();
835 }
836 
837 
838 //----------------------------------------------------------------------//
839 
840 template <typename MatrixType, typename PreconditionerType, typename VectorType>
841 inline
843 (const unsigned int steps,
844  const bool variable,
845  const bool symmetric,
846  const bool transpose)
847  :
848  MGSmoother<VectorType>(steps, variable, symmetric, transpose)
849 {}
850 
851 
852 
853 template <typename MatrixType, typename PreconditionerType, typename VectorType>
854 inline void
856 {
857  smoothers.clear();
858 
859  unsigned int i=matrices.min_level(),
860  max_level=matrices.max_level();
861  for (; i<=max_level; ++i)
862  matrices[i]=0;
863 }
864 
865 
866 
867 template <typename MatrixType, typename PreconditionerType, typename VectorType>
868 template <typename MatrixType2>
869 inline void
872  const typename PreconditionerType::AdditionalData &data)
873 {
874  const unsigned int min = m.min_level();
875  const unsigned int max = m.max_level();
876 
877  matrices.resize(min, max);
878  smoothers.resize(min, max);
879 
880  for (unsigned int i=min; i<=max; ++i)
881  {
882  matrices[i] = &m[i];
883  smoothers[i].initialize(m[i], data);
884  }
885 }
886 
887 
888 
889 template <typename MatrixType, typename PreconditionerType, typename VectorType>
890 template <typename MatrixType2, class DATA>
891 inline void
894  const MGLevelObject<DATA> &data)
895 {
896  const unsigned int min = m.min_level();
897  const unsigned int max = m.max_level();
898 
899  Assert (data.min_level() == min,
900  ExcDimensionMismatch(data.min_level(), min));
901  Assert (data.max_level() == max,
902  ExcDimensionMismatch(data.max_level(), max));
903 
904  matrices.resize(min, max);
905  smoothers.resize(min, max);
906 
907  for (unsigned int i=min; i<=max; ++i)
908  {
909  matrices[i] = &m[i];
910  smoothers[i].initialize(m[i], data[i]);
911  }
912 }
913 
914 
915 
916 template <typename MatrixType, typename PreconditionerType, typename VectorType>
917 template <typename MatrixType2, class DATA>
918 inline void
921  const DATA &data,
922  const unsigned int row,
923  const unsigned int col)
924 {
925  const unsigned int min = m.min_level();
926  const unsigned int max = m.max_level();
927 
928  matrices.resize(min, max);
929  smoothers.resize(min, max);
930 
931  for (unsigned int i=min; i<=max; ++i)
932  {
933  matrices[i] = &(m[i].block(row, col));
934  smoothers[i].initialize(m[i].block(row, col), data);
935  }
936 }
937 
938 
939 
940 template <typename MatrixType, typename PreconditionerType, typename VectorType>
941 template <typename MatrixType2, class DATA>
942 inline void
945  const MGLevelObject<DATA> &data,
946  const unsigned int row,
947  const unsigned int col)
948 {
949  const unsigned int min = m.min_level();
950  const unsigned int max = m.max_level();
951 
952  Assert (data.min_level() == min,
953  ExcDimensionMismatch(data.min_level(), min));
954  Assert (data.max_level() == max,
955  ExcDimensionMismatch(data.max_level(), max));
956 
957  matrices.resize(min, max);
958  smoothers.resize(min, max);
959 
960  for (unsigned int i=min; i<=max; ++i)
961  {
962  matrices[i] = &(m[i].block(row, col));
963  smoothers[i].initialize(m[i].block(row, col), data[i]);
964  }
965 }
966 
967 
968 
969 template <typename MatrixType, typename PreconditionerType, typename VectorType>
970 inline void
972 (const unsigned int level,
973  VectorType &u,
974  const VectorType &rhs) const
975 {
976  unsigned int maxlevel = matrices.max_level();
977  unsigned int steps2 = this->steps;
978 
979  if (this->variable)
980  steps2 *= (1<<(maxlevel-level));
981 
982  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
983  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
984 
985  r->reinit(u,true);
986  d->reinit(u,true);
987 
988  bool T = this->transpose;
989  if (this->symmetric && (steps2 % 2 == 0))
990  T = false;
991  if (this->debug > 0)
992  deallog << 'S' << level << ' ';
993 
994  for (unsigned int i=0; i<steps2; ++i)
995  {
996  if (T)
997  {
998  if (this->debug > 0)
999  deallog << 'T';
1000  if (i == 0 && u.all_zero())
1001  *r = rhs;
1002  else
1003  {
1004  matrices[level].Tvmult(*r,u);
1005  r->sadd(-1.,1.,rhs);
1006  }
1007  if (this->debug > 2)
1008  deallog << ' ' << r->l2_norm() << ' ';
1009  smoothers[level].Tvmult(*d, *r);
1010  if (this->debug > 1)
1011  deallog << ' ' << d->l2_norm() << ' ';
1012  }
1013  else
1014  {
1015  if (this->debug > 0)
1016  deallog << 'N';
1017  if (i == 0 && u.all_zero())
1018  *r = rhs;
1019  else
1020  {
1021  matrices[level].vmult(*r,u);
1022  r->sadd(-1.,rhs);
1023  }
1024  if (this->debug > 2)
1025  deallog << ' ' << r->l2_norm() << ' ';
1026  smoothers[level].vmult(*d, *r);
1027  if (this->debug > 1)
1028  deallog << ' ' << d->l2_norm() << ' ';
1029  }
1030  u += *d;
1031  if (this->symmetric)
1032  T = !T;
1033  }
1034  if (this->debug > 0)
1035  deallog << std::endl;
1036 }
1037 
1038 
1039 
1040 template <typename MatrixType, typename PreconditionerType, typename VectorType>
1041 inline
1042 std::size_t
1044 memory_consumption () const
1045 {
1046  return sizeof(*this)
1047  + matrices.memory_consumption()
1048  + smoothers.memory_consumption()
1049  + this->vector_memory.memory_consumption();
1050 }
1051 
1052 
1053 #endif // DOXYGEN
1054 
1055 DEAL_II_NAMESPACE_CLOSE
1056 
1057 #endif
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void set_transpose(const bool)
bool symmetric
Definition: mg_smoother.h:109
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
void set_symmetric(const bool)
Definition: mg.h:81
std::size_t memory_consumption() const
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
virtual void clear()
void set_variable(const bool)
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
MGLevelObject< RelaxationType > smoothers
Definition: mg_smoother.h:347
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
std::size_t memory_consumption() const
void set_debug(const unsigned int level)
MGLevelObject< PreconditionerType > smoothers
Definition: mg_smoother.h:480
VectorizedArray< Number > min(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
#define Assert(cond, exc)
Definition: exceptions.h:294
unsigned int max_level() const
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< PointerMatrix< MatrixType, VectorType > > matrices
Definition: mg_smoother.h:359
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
GrowingVectorMemory< VectorType > vector_memory
Definition: mg_smoother.h:91
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
void set_steps(const unsigned int)
std::size_t memory_consumption() const
MGLevelObject< PointerMatrix< MatrixType, VectorType > > matrices
Definition: mg_smoother.h:492
std::size_t memory_consumption() const
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
unsigned int debug
Definition: mg_smoother.h:120
bool transpose
Definition: mg_smoother.h:115
unsigned int min_level() const
bool variable
Definition: mg_smoother.h:103
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
VectorizedArray< Number > max(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
unsigned int steps
Definition: mg_smoother.h:97
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())