dune-pdelab  2.5-dev
ovlpistlsolverbackend.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
5 
6 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
7 #include "ensureistlinclude.hh"
8 
9 #include <dune/common/deprecated.hh>
10 #include <dune/common/parallel/mpihelper.hh>
11 
12 #include <dune/istl/owneroverlapcopy.hh>
13 #include <dune/istl/solvercategory.hh>
14 #include <dune/istl/operators.hh>
15 #include <dune/istl/solvers.hh>
16 #include <dune/istl/preconditioners.hh>
17 #include <dune/istl/scalarproducts.hh>
18 #include <dune/istl/paamg/amg.hh>
19 #include <dune/istl/paamg/pinfo.hh>
20 #include <dune/istl/io.hh>
21 #include <dune/istl/superlu.hh>
22 
29 
30 namespace Dune {
31  namespace PDELab {
32 
36 
37  //========================================================
38  // Generic support for overlapping grids
39  // (need to be used with appropriate constraints)
40  //========================================================
41 
42  // operator that resets result to zero at constrained DOFS
43  template<class CC, class M, class X, class Y>
45  : public Dune::AssembledLinearOperator<M,X,Y>
46  {
47  public:
49  typedef M matrix_type;
50  typedef X domain_type;
51  typedef Y range_type;
52  typedef typename X::ElementType field_type;
53 
54  //redefine the category, that is the only difference
55  enum {category=Dune::SolverCategory::overlapping};
56 
57  OverlappingOperator (const CC& cc_, const M& A)
58  : cc(cc_), _A_(A)
59  {}
60 
62  virtual void apply (const domain_type& x, range_type& y) const
63  {
64  using Backend::native;
65  native(_A_).mv(native(x),native(y));
67  }
68 
70  virtual void applyscaleadd (field_type alpha, const domain_type& x, range_type& y) const
71  {
72  using Backend::native;
73  native(_A_).usmv(alpha,native(x),native(y));
75  }
76 
78  virtual const M& getmat () const
79  {
80  return _A_;
81  }
82 
83  private:
84  const CC& cc;
85  const M& _A_;
86  };
87 
88  // new scalar product assuming at least overlap 1
89  // uses unique partitioning of nodes for parallelization
90  template<class GFS, class X>
92  : public Dune::ScalarProduct<X>
93  {
94  public:
96  typedef X domain_type;
97  typedef typename X::ElementType field_type;
98 
100  enum {category=Dune::SolverCategory::overlapping};
101 
104  OverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
105  : gfs(gfs_), helper(helper_)
106  {}
107 
108 
113  virtual field_type dot (const X& x, const X& y)
114  {
115  // do local scalar product on unique partition
116  field_type sum = helper.disjointDot(x,y);
117 
118  // do global communication
119  return gfs.gridView().comm().sum(sum);
120  }
121 
125  virtual double norm (const X& x)
126  {
127  return sqrt(static_cast<double>(this->dot(x,x)));
128  }
129 
130  private:
131  const GFS& gfs;
132  const ISTL::ParallelHelper<GFS>& helper;
133  };
134 
135  // wrapped sequential preconditioner
136  template<class CC, class GFS, class P>
138  : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<GFS,typename P::domain_type::field_type>,
139  Dune::PDELab::Backend::Vector<GFS,typename P::range_type::field_type>>
140  {
141  public:
146 
147  // define the category
148  enum {
150  category=Dune::SolverCategory::overlapping
151  };
152 
154  OverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
155  const ISTL::ParallelHelper<GFS>& helper_)
156  : gfs(gfs_), prec(prec_), cc(cc_), helper(helper_)
157  {}
158 
162  virtual void pre (domain_type& x, range_type& b)
163  {
164  prec.pre(x,b);
165  }
166 
170  virtual void apply (domain_type& v, const range_type& d)
171  {
172  range_type dd(d);
173  set_constrained_dofs(cc,0.0,dd);
174  prec.apply(Backend::native(v),Backend::native(dd));
176  if (gfs.gridView().comm().size()>1)
177  gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
178  }
179 
183  virtual void post (domain_type& x)
184  {
185  prec.post(Backend::native(x));
186  }
187 
188  private:
189  const GFS& gfs;
190  P& prec;
191  const CC& cc;
192  const ISTL::ParallelHelper<GFS>& helper;
193  };
194 
195 
196 #if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
197  // exact subdomain solves with UMFPack as preconditioner
198  template<class GFS, class M, class X, class Y>
199  class UMFPackSubdomainSolver : public Dune::Preconditioner<X,Y>
200  {
201  typedef Backend::Native<M> ISTLM;
202 
203  public:
205  typedef X domain_type;
207  typedef Y range_type;
209  typedef typename X::ElementType field_type;
210 
211 
212  // define the category
213  enum {
215  category=Dune::SolverCategory::overlapping
216  };
217 
224  UMFPackSubdomainSolver (const GFS& gfs_, const M& A_)
225  : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
226  {}
227 
231  virtual void pre (X& x, Y& b) {}
232 
236  virtual void apply (X& v, const Y& d)
237  {
238  Dune::InverseOperatorResult stat;
239  Y b(d); // need copy, since solver overwrites right hand side
240  solver.apply(Backend::native(v),Backend::native(b),stat);
241  if (gfs.gridView().comm().size()>1)
242  {
243  AddDataHandle<GFS,X> adddh(gfs,v);
244  gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
245  }
246  }
247 
251  virtual void post (X& x) {}
252 
253  private:
254  const GFS& gfs;
255  Dune::UMFPack<ISTLM> solver;
256  };
257 #endif
258 
259 #if HAVE_SUPERLU
260  // exact subdomain solves with SuperLU as preconditioner
261  template<class GFS, class M, class X, class Y>
262  class SuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
263  {
264  typedef Backend::Native<M> ISTLM;
265 
266  public:
268  typedef X domain_type;
270  typedef Y range_type;
272  typedef typename X::ElementType field_type;
273 
274 
275  // define the category
276  enum {
278  category=Dune::SolverCategory::overlapping
279  };
280 
287  SuperLUSubdomainSolver (const GFS& gfs_, const M& A_)
288  : gfs(gfs_), solver(Backend::native(A_),false) // this does the decomposition
289  {}
290 
294  virtual void pre (X& x, Y& b) {}
295 
299  virtual void apply (X& v, const Y& d)
300  {
301  Dune::InverseOperatorResult stat;
302  Y b(d); // need copy, since solver overwrites right hand side
303  solver.apply(Backend::native(v),Backend::native(b),stat);
304  if (gfs.gridView().comm().size()>1)
305  {
306  AddDataHandle<GFS,X> adddh(gfs,v);
307  gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
308  }
309  }
310 
314  virtual void post (X& x) {}
315 
316  private:
317  const GFS& gfs;
318  Dune::SuperLU<ISTLM> solver;
319  };
320 
321  // exact subdomain solves with SuperLU as preconditioner
322  template<class GFS, class M, class X, class Y>
323  class RestrictedSuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
324  {
325  typedef typename M::Container ISTLM;
326 
327  public:
329  typedef X domain_type;
331  typedef Y range_type;
333  typedef typename X::ElementType field_type;
334 
335 
336  // define the category
337  enum {
339  category=Dune::SolverCategory::overlapping
340  };
341 
349  RestrictedSuperLUSubdomainSolver (const GFS& gfs_, const M& A_,
350  const ISTL::ParallelHelper<GFS>& helper_)
351  : gfs(gfs_), solver(Backend::native(A_),false), helper(helper_) // this does the decomposition
352  {}
353 
357  virtual void pre (X& x, Y& b) {}
358 
362  virtual void apply (X& v, const Y& d)
363  {
364  using Backend::native;
365  Dune::InverseOperatorResult stat;
366  Y b(d); // need copy, since solver overwrites right hand side
367  solver.apply(native(v),native(b),stat);
368  if (gfs.gridView().comm().size()>1)
369  {
370  helper.maskForeignDOFs(native(v));
371  AddDataHandle<GFS,X> adddh(gfs,v);
372  gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
373  }
374  }
375 
379  virtual void post (X& x) {}
380 
381  private:
382  const GFS& gfs;
383  Dune::SuperLU<ISTLM> solver;
384  const ISTL::ParallelHelper<GFS>& helper;
385  };
386 #endif
387 
388  template<typename GFS>
390  {
391  public:
393  : gfs(gfs_), helper(gfs_)
394  {}
395 
400  template<typename X>
401  typename X::ElementType dot (const X& x, const X& y) const
402  {
403  // do local scalar product on unique partition
404  typename X::ElementType sum = helper.disjointDot(x,y);
405 
406  // do global communication
407  return gfs.gridView().comm().sum(sum);
408  }
409 
413  template<typename X>
414  typename Dune::template FieldTraits<typename X::ElementType >::real_type norm (const X& x) const
415  {
416  using namespace std;
417  return sqrt(static_cast<double>(this->dot(x,x)));
418  }
419 
421  {
422  return helper;
423  }
424 
425  // need also non-const version;
426  ISTL::ParallelHelper<GFS>& parallelHelper() // P.B.: needed for createIndexSetAndProjectForAMG
427  {
428  return helper;
429  }
430 
431  private:
432  const GFS& gfs;
434  };
435 
436 
437  template<typename GFS, typename X>
439  : public ScalarProduct<X>
440  {
441  public:
442  enum {category=Dune::SolverCategory::overlapping};
444  : implementation(implementation_)
445  {}
446 
447  virtual typename X::Container::field_type dot(const X& x, const X& y)
448  {
449  return implementation.dot(x,y);
450  }
451 
452  virtual typename X::Container::field_type norm (const X& x)
453  {
454  using namespace std;
455  return sqrt(static_cast<double>(this->dot(x,x)));
456  }
457 
458  private:
459  const OVLPScalarProductImplementation<GFS>& implementation;
460  };
461 
462  template<class GFS, class C,
463  template<class,class,class,int> class Preconditioner,
464  template<class> class Solver>
467  {
468  public:
477  ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
478  int steps_=5, int verbose_=1)
479  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
480  {}
481 
489  template<class M, class V, class W>
490  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
491  {
492  using Backend::Native;
493  using Backend::native;
494  typedef OverlappingOperator<C,M,V,W> POP;
495  POP pop(c,A);
496  typedef OVLPScalarProduct<GFS,V> PSP;
497  PSP psp(*this);
498  typedef Preconditioner<
499  Native<M>,
500  Native<V>,
501  Native<W>,
502  1
503  > SeqPrec;
504  SeqPrec seqprec(native(A),steps,1.0);
506  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
507  int verb=0;
508  if (gfs.gridView().comm().rank()==0) verb=verbose;
509  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
510  Dune::InverseOperatorResult stat;
511  solver.apply(z,r,stat);
512  res.converged = stat.converged;
513  res.iterations = stat.iterations;
514  res.elapsed = stat.elapsed;
515  res.reduction = stat.reduction;
516  res.conv_rate = stat.conv_rate;
517  }
518  private:
519  const GFS& gfs;
520  const C& c;
521  unsigned maxiter;
522  int steps;
523  int verbose;
524  };
525 
526  // Base class for ILU0 as preconditioner
527  template<class GFS, class C,
528  template<class> class Solver>
531  {
532  public:
540  ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
541  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
542  {}
543 
551  template<class M, class V, class W>
552  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
553  {
554  using Backend::Native;
555  using Backend::native;
556  typedef OverlappingOperator<C,M,V,W> POP;
557  POP pop(c,A);
558  typedef OVLPScalarProduct<GFS,V> PSP;
559  PSP psp(*this);
560  typedef SeqILU0<
561  Native<M>,
562  Native<V>,
563  Native<W>,
564  1
565  > SeqPrec;
566  SeqPrec seqprec(native(A),1.0);
568  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
569  int verb=0;
570  if (gfs.gridView().comm().rank()==0) verb=verbose;
571  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
572  Dune::InverseOperatorResult stat;
573  solver.apply(z,r,stat);
574  res.converged = stat.converged;
575  res.iterations = stat.iterations;
576  res.elapsed = stat.elapsed;
577  res.reduction = stat.reduction;
578  res.conv_rate = stat.conv_rate;
579  }
580  private:
581  const GFS& gfs;
582  const C& c;
583  unsigned maxiter;
584  int steps;
585  int verbose;
586  };
587 
588  // Base class for ILUn as preconditioner
589  template<class GFS, class C,
590  template<class> class Solver>
593  {
594  public:
603  ISTLBackend_OVLP_ILUn_Base (const GFS& gfs_, const C& c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
604  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), n(n_), maxiter(maxiter_), verbose(verbose_)
605  {}
606 
614  template<class M, class V, class W>
615  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
616  {
617  using Backend::Native;
618  using Backend::native;
619  typedef OverlappingOperator<C,M,V,W> POP;
620  POP pop(c,A);
621  typedef OVLPScalarProduct<GFS,V> PSP;
622  PSP psp(*this);
623  typedef SeqILUn<
624  Native<M>,
625  Native<V>,
626  Native<W>,
627  1
628  > SeqPrec;
629  SeqPrec seqprec(native(A),n,1.0);
631  WPREC wprec(gfs,seqprec,c,this->parallelHelper());
632  int verb=0;
633  if (gfs.gridView().comm().rank()==0) verb=verbose;
634  Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
635  Dune::InverseOperatorResult stat;
636  solver.apply(z,r,stat);
637  res.converged = stat.converged;
638  res.iterations = stat.iterations;
639  res.elapsed = stat.elapsed;
640  res.reduction = stat.reduction;
641  res.conv_rate = stat.conv_rate;
642  }
643  private:
644  const GFS& gfs;
645  const C& c;
646  int n;
647  unsigned maxiter;
648  int steps;
649  int verbose;
650  };
651 
654 
660  template<class GFS, class CC>
662  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
663  {
664  public:
673  ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
674  int steps=5, int verbose=1)
675  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
676  {}
677  };
683  template<class GFS, class CC>
685  : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
686  {
687  public:
695  ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
696  : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
697  {}
698  };
704  template<class GFS, class CC>
706  : public ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>
707  {
708  public:
717  ISTLBackend_OVLP_BCGS_ILUn (const GFS& gfs, const CC& cc, int n=1, unsigned maxiter=5000, int verbose=1)
718  : ISTLBackend_OVLP_ILUn_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, n, maxiter, verbose)
719  {}
720  };
726  template<class GFS, class CC>
728  : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
729  {
730  public:
739  ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
740  int steps=5, int verbose=1)
741  : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
742  {}
743  };
744 
750  template<class GFS, class CC>
753  {
754  public:
762  ISTLBackend_OVLP_GMRES_ILU0 (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000, int verbose_=1,
763  int restart_ = 20)
764  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), cc(cc_), maxiter(maxiter_), verbose(verbose_),
765  restart(restart_)
766  {}
767 
774  template<class M, class V, class W>
775  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
776  {
777  using Backend::Native;
778  using Backend::native;
779  typedef OverlappingOperator<CC,M,V,W> POP;
780  POP pop(cc,A);
781  typedef OVLPScalarProduct<GFS,V> PSP;
782  PSP psp(*this);
783  typedef SeqILU0<
784  Native<M>,
785  Native<V>,
786  Native<W>,
787  1
788  > SeqPrec;
789  SeqPrec seqprec(native(A),1.0);
791  WPREC wprec(gfs,seqprec,cc,this->parallelHelper());
792  int verb=0;
793  if (gfs.gridView().comm().rank()==0) verb=verbose;
794  RestartedGMResSolver<V> solver(pop,psp,wprec,reduction,restart,maxiter,verb);
795  Dune::InverseOperatorResult stat;
796  solver.apply(z,r,stat);
797  res.converged = stat.converged;
798  res.iterations = stat.iterations;
799  res.elapsed = stat.elapsed;
800  res.reduction = stat.reduction;
801  res.conv_rate = stat.conv_rate;
802  }
803 
804  private:
805  const GFS& gfs;
806  const CC& cc;
807  unsigned maxiter;
808  int steps;
809  int verbose;
810  int restart;
811  };
812 
814 
815  template<class GFS, class C, template<typename> class Solver>
818  {
819  public:
827  ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
828  int verbose_=1)
829  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
830  {}
831 
839  template<class M, class V, class W>
840  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
841  {
842  typedef OverlappingOperator<C,M,V,W> POP;
843  POP pop(c,A);
844  typedef OVLPScalarProduct<GFS,V> PSP;
845  PSP psp(*this);
846 #if HAVE_SUPERLU
847  typedef SuperLUSubdomainSolver<GFS,M,V,W> PREC;
848  PREC prec(gfs,A);
849  int verb=0;
850  if (gfs.gridView().comm().rank()==0) verb=verbose;
851  Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
852  Dune::InverseOperatorResult stat;
853  solver.apply(z,r,stat);
854  res.converged = stat.converged;
855  res.iterations = stat.iterations;
856  res.elapsed = stat.elapsed;
857  res.reduction = stat.reduction;
858  res.conv_rate = stat.conv_rate;
859 #else
860  std::cout << "No superLU support, please install and configure it." << std::endl;
861 #endif
862  }
863 
864  private:
865  const GFS& gfs;
866  const C& c;
867  unsigned maxiter;
868  int verbose;
869  };
870 
872 
873  template<class GFS, class C, template<typename> class Solver>
876  {
877  public:
885  ISTLBackend_OVLP_UMFPack_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
886  int verbose_=1)
887  : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
888  {}
889 
897  template<class M, class V, class W>
898  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
899  {
900  typedef OverlappingOperator<C,M,V,W> POP;
901  POP pop(c,A);
902  typedef OVLPScalarProduct<GFS,V> PSP;
903  PSP psp(*this);
904 #if HAVE_SUITESPARSE_UMFPACK || DOXYGEN
906  PREC prec(gfs,A);
907  int verb=0;
908  if (gfs.gridView().comm().rank()==0) verb=verbose;
909  Solver<V> solver(pop,psp,prec,reduction,maxiter,verb);
910  Dune::InverseOperatorResult stat;
911  solver.apply(z,r,stat);
912  res.converged = stat.converged;
913  res.iterations = stat.iterations;
914  res.elapsed = stat.elapsed;
915  res.reduction = stat.reduction;
916  res.conv_rate = stat.conv_rate;
917 #else
918  std::cout << "No UMFPack support, please install and configure it." << std::endl;
919 #endif
920  }
921 
922  private:
923  const GFS& gfs;
924  const C& c;
925  unsigned maxiter;
926  int verbose;
927  };
928 
931 
936  template<class GFS, class CC>
938  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
939  {
940  public:
941 
949  ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
950  int verbose_=1)
951  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
952  {}
953  };
954 
960  template<class GFS, class CC>
962  : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
963  {
964  public:
965 
973  ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
974  unsigned maxiter_=5000,
975  int verbose_=1)
976  : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
977  {}
978  };
979 
985  template<class GFS, class CC>
987  : public ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>
988  {
989  public:
990 
998  ISTLBackend_OVLP_CG_UMFPack (const GFS& gfs_, const CC& cc_,
999  unsigned maxiter_=5000,
1000  int verbose_=1)
1001  : ISTLBackend_OVLP_UMFPack_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
1002  {}
1003  };
1004 
1005 
1009  template<class GFS>
1011  : public LinearResultStorage
1012  {
1013  public:
1018  explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
1019  : gfs(gfs_)
1020  {}
1021 
1023  : gfs(other_.gfs)
1024  {}
1025 
1030  template<class V>
1031  typename V::ElementType norm(const V& v) const
1032  {
1033  static_assert
1035  "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
1036  "neccessary, so we skipped the implementation. If you have a "
1037  "scenario where you need it, please implement it or report back to "
1038  "us.");
1039  }
1040 
1048  template<class M, class V, class W>
1049  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
1050  {
1051  using Backend::Native;
1052  using Backend::native;
1053  Dune::SeqJac<
1054  Native<M>,
1055  Native<V>,
1056  Native<W>
1057  > jac(native(A),1,1.0);
1058  jac.pre(native(z),native(r));
1059  jac.apply(native(z),native(r));
1060  jac.post(native(z));
1061  if (gfs.gridView().comm().size()>1)
1062  {
1063  CopyDataHandle<GFS,V> copydh(gfs,z);
1064  gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
1065  }
1066  res.converged = true;
1067  res.iterations = 1;
1068  res.elapsed = 0.0;
1069  res.reduction = static_cast<double>(reduction);
1070  res.conv_rate = static_cast<double>(reduction); // pow(reduction,1.0/1)
1071  }
1072 
1073  private:
1074  const GFS& gfs;
1075  };
1077 
1078  template<class GO, int s, template<class,class,class,int> class Preconditioner,
1079  template<class> class Solver>
1081  {
1082  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1084  typedef typename GO::Traits::Jacobian M;
1085  typedef Backend::Native<M> MatrixType;
1086  typedef typename GO::Traits::Domain V;
1087  typedef Backend::Native<V> VectorType;
1089 #if HAVE_MPI
1090  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
1091  typedef Dune::BlockPreconditioner<VectorType,VectorType,Comm,Smoother> ParSmoother;
1092  typedef Dune::OverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
1093 #else
1094  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
1095  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
1096 #endif
1097  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
1098  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
1099 
1100  typedef typename V::ElementType RF;
1101 
1102  public:
1103 
1107  typedef Dune::Amg::Parameters Parameters;
1108 
1109  public:
1110  ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
1111  int verbose_=1, bool reuse_=false,
1112  bool usesuperlu_=true)
1113  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000),
1114  verbose(verbose_), reuse(reuse_), firstapply(true),
1115  usesuperlu(usesuperlu_)
1116  {
1117  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
1118  params.setDebugLevel(verbose_);
1119 #if !HAVE_SUPERLU
1120  if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
1121  {
1122  std::cout << "WARNING: You are using AMG without SuperLU!"
1123  << " Please consider installing SuperLU,"
1124  << " or set the usesuperlu flag to false"
1125  << " to suppress this warning." << std::endl;
1126  }
1127 #endif
1128  }
1129 
1134  void setParameters(const Parameters& params_)
1135  {
1136  params = params_;
1137  }
1138 
1146  const Parameters& parameters() const
1147  {
1148  return params;
1149  }
1150 
1152  void setReuse(bool reuse_)
1153  {
1154  reuse = reuse_;
1155  }
1156 
1158  bool getReuse() const
1159  {
1160  return reuse;
1161  }
1162 
1167  typename V::ElementType norm (const V& v) const
1168  {
1169  typedef OverlappingScalarProduct<GFS,V> PSP;
1170  PSP psp(gfs,phelper);
1171  return psp.norm(v);
1172  }
1173 
1181  void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
1182  {
1183  Timer watch;
1184  Comm oocc(gfs.gridView().comm());
1185  MatrixType& mat=Backend::native(A);
1186  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
1187  Dune::Amg::FirstDiagonal> > Criterion;
1188 #if HAVE_MPI
1189  phelper.createIndexSetAndProjectForAMG(A, oocc);
1190  Operator oop(mat, oocc);
1191  Dune::OverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
1192 #else
1193  Operator oop(mat);
1194  Dune::SeqScalarProduct<VectorType> sp;
1195 #endif
1196  SmootherArgs smootherArgs;
1197  smootherArgs.iterations = 1;
1198  smootherArgs.relaxationFactor = 1;
1199  Criterion criterion(params);
1200  stats.tprepare=watch.elapsed();
1201  watch.reset();
1202 
1203  int verb=0;
1204  if (gfs.gridView().comm().rank()==0) verb=verbose;
1205  //only construct a new AMG if the matrix changes
1206  if (reuse==false || firstapply==true){
1207  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1208  firstapply = false;
1209  stats.tsetup = watch.elapsed();
1210  stats.levels = amg->maxlevels();
1211  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1212  }
1213  watch.reset();
1214  Solver<VectorType> solver(oop,sp,*amg,RF(reduction),maxiter,verb);
1215  Dune::InverseOperatorResult stat;
1216 
1217  solver.apply(Backend::native(z),Backend::native(r),stat);
1218  stats.tsolve= watch.elapsed();
1219  res.converged = stat.converged;
1220  res.iterations = stat.iterations;
1221  res.elapsed = stat.elapsed;
1222  res.reduction = stat.reduction;
1223  res.conv_rate = stat.conv_rate;
1224  }
1225 
1231  {
1232  return stats;
1233  }
1234 
1235  private:
1236  const GFS& gfs;
1237  PHELPER phelper;
1238  unsigned maxiter;
1239  Parameters params;
1240  int verbose;
1241  bool reuse;
1242  bool firstapply;
1243  bool usesuperlu;
1244  shared_ptr<AMG> amg;
1245  ISTLAMGStatistics stats;
1246  };
1247 
1250 
1257  template<class GO, int s=96>
1259  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1260  {
1261  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1262  public:
1272  ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1273  int verbose_=1, bool reuse_=false,
1274  bool usesuperlu_=true)
1275  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1276  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1277  {}
1278  };
1279 
1286  template<class GO, int s=96>
1288  : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1289  {
1290  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1291  public:
1301  ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
1302  int verbose_=1, bool reuse_=false,
1303  bool usesuperlu_=true)
1304  : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1305  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1306  {}
1307  };
1308 
1315  template<class GO, int s=96>
1317  : public ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1318  {
1319  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
1320  public:
1330  ISTLBackend_BCGS_AMG_ILU0(const GFS& gfs_, unsigned maxiter_=5000,
1331  int verbose_=1, bool reuse_=false,
1332  bool usesuperlu_=true)
1333  : ISTLBackend_AMG<GO, s, Dune::SeqILU0, Dune::BiCGSTABSolver>
1334  (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
1335  {}
1336  };
1337 
1339 
1341 
1342  } // namespace PDELab
1343 } // namespace Dune
1344 
1345 #endif // DUNE_PDELAB_BACKEND_ISTL_OVLPISTLSOLVERBACKEND_HH
M matrix_type
export types
Definition: ovlpistlsolverbackend.hh:49
ISTLBackend_OVLP_ExplicitDiagonal(const ISTLBackend_OVLP_ExplicitDiagonal &other_)
Definition: ovlpistlsolverbackend.hh:1022
Overlapping parallel CG solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:961
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
virtual const M & getmat() const
get matrix via *
Definition: ovlpistlsolverbackend.hh:78
OverlappingOperator(const CC &cc_, const M &A)
Definition: ovlpistlsolverbackend.hh:57
virtual void post(domain_type &x)
Clean up.
Definition: ovlpistlsolverbackend.hh:183
virtual void post(X &x)
Clean up.
Definition: ovlpistlsolverbackend.hh:251
OVLPScalarProduct(const OVLPScalarProductImplementation< GFS > &implementation_)
Definition: ovlpistlsolverbackend.hh:443
ISTLBackend_BCGS_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1301
void setParameters(const Parameters &params_)
set AMG parameters
Definition: ovlpistlsolverbackend.hh:1134
ISTLBackend_OVLP_ILUn_Base(const GFS &gfs_, const C &c_, int n_=1, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:603
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper...
Definition: backend/interface.hh:176
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: ovlpistlsolverbackend.hh:236
Class providing some statistics of the AMG solver.
Definition: seqistlsolverbackend.hh:544
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: ovlpistlsolverbackend.hh:1287
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: ovlpistlsolverbackend.hh:1230
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: ovlpistlsolverbackend.hh:1010
ISTL::ParallelHelper< GFS > & parallelHelper()
Definition: ovlpistlsolverbackend.hh:426
Definition: ovlpistlsolverbackend.hh:55
ISTLBackend_OVLP_UMFPack_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:885
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:1167
ISTLBackend_OVLP_BCGS_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:949
Dune::PDELab::Backend::Vector< GFS, typename P::domain_type::field_type > domain_type
The domain type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:143
OVLPScalarProductImplementation(const GFS &gfs_)
Definition: ovlpistlsolverbackend.hh:392
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: ovlpistlsolverbackend.hh:1158
ISTLBackend_OVLP_BCGS_SSORk(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int steps=5, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:673
Definition: ovlpistlsolverbackend.hh:137
Definition: solver.hh:42
Definition: parallelhelper.hh:53
const std::string s
Definition: function.hh:1101
ISTLBackend_OVLP_BCGS_ILUn(const GFS &gfs, const CC &cc, int n=1, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:717
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: ovlpistlsolverbackend.hh:1146
ISTLBackend_OVLP_ILU0_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:540
virtual void applyscaleadd(field_type alpha, const domain_type &x, range_type &y) const
apply operator to x, scale and add:
Definition: ovlpistlsolverbackend.hh:70
virtual void pre(domain_type &x, range_type &b)
Prepare the preconditioner.
Definition: ovlpistlsolverbackend.hh:162
ISTLBackend_OVLP_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:477
virtual double norm(const X &x)
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: ovlpistlsolverbackend.hh:125
OverlappingWrappedPreconditioner(const GFS &gfs_, P &prec_, const CC &cc_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor.
Definition: ovlpistlsolverbackend.hh:154
X domain_type
Definition: ovlpistlsolverbackend.hh:50
Y range_type
The range type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:207
Definition: ovlpistlsolverbackend.hh:529
X::ElementType dot(const X &x, const X &y) const
Dot product of two vectors. It is assumed that the vectors are consistent on the interior+border part...
Definition: ovlpistlsolverbackend.hh:401
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Definition: ovlpistlsolverbackend.hh:389
Definition: genericdatahandle.hh:622
ISTLBackend_OVLP_GMRES_ILU0(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1, int restart_=20)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:762
ISTLBackend_OVLP_SuperLU_Base(const GFS &gfs_, const C &c_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:827
ISTLBackend_OVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:1018
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: ovlpistlsolverbackend.hh:231
OverlappingScalarProduct(const GFS &gfs_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: ovlpistlsolverbackend.hh:104
Overlapping parallel CGS solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:727
X domain_type
export types
Definition: ovlpistlsolverbackend.hh:96
Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR...
Definition: ovlpistlsolverbackend.hh:1258
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:840
Overlapping parallel BiCGStab solver with SSOR preconditioner.
Definition: ovlpistlsolverbackend.hh:661
virtual void apply(const domain_type &x, range_type &y) const
apply operator to x:
Definition: ovlpistlsolverbackend.hh:62
Overlapping parallel CG solver with UMFPack preconditioner.
Definition: ovlpistlsolverbackend.hh:986
virtual X::Container::field_type dot(const X &x, const X &y)
Definition: ovlpistlsolverbackend.hh:447
ISTLBackend_OVLP_CG_UMFPack(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:998
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
virtual X::Container::field_type norm(const X &x)
Definition: ovlpistlsolverbackend.hh:452
X::ElementType field_type
Definition: ovlpistlsolverbackend.hh:52
Definition: ovlpistlsolverbackend.hh:91
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:552
Definition: ovlpistlsolverbackend.hh:816
Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by ILU0.
Definition: ovlpistlsolverbackend.hh:1316
ISTLBackend_OVLP_BCGS_ILU0(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:695
Definition: ovlpistlsolverbackend.hh:1080
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:775
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:1181
virtual field_type dot(const X &x, const X &y)
Dot product of two vectors. It is assumed that the vectors are consistent on the interior+border part...
Definition: ovlpistlsolverbackend.hh:113
Definition: ovlpistlsolverbackend.hh:465
Definition: ovlpistlsolverbackend.hh:199
Dune::template FieldTraits< typename X::ElementType >::real_type norm(const X &x) const
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: ovlpistlsolverbackend.hh:414
Definition: ovlpistlsolverbackend.hh:44
X::ElementType field_type
Definition: ovlpistlsolverbackend.hh:97
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: ovlpistlsolverbackend.hh:1152
STL namespace.
ISTLBackend_OVLP_CG_SuperLU(const GFS &gfs_, const CC &cc_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:973
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename W::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:1049
UMFPackSubdomainSolver(const GFS &gfs_, const M &A_)
Constructor.
Definition: ovlpistlsolverbackend.hh:224
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:705
Dune::Amg::Parameters Parameters
Parameters object to customize matrix hierachy building.
Definition: ovlpistlsolverbackend.hh:1107
X::ElementType field_type
The field type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:209
Definition: ovlpistlsolverbackend.hh:591
const ISTL::ParallelHelper< GFS > & parallelHelper() const
Definition: ovlpistlsolverbackend.hh:420
ISTLBackend_OVLP_CG_SSORk(const GFS &gfs, const CC &cc, unsigned maxiter=5000, int steps=5, int verbose=1)
make a linear solver object
Definition: ovlpistlsolverbackend.hh:739
ISTLBackend_AMG(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: ovlpistlsolverbackend.hh:1110
Definition: genericdatahandle.hh:685
Definition: ovlpistlsolverbackend.hh:438
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:615
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:798
Overlapping parallel BiCGStab solver with SuperLU preconditioner.
Definition: ovlpistlsolverbackend.hh:937
ISTLBackend_BCGS_AMG_ILU0(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1330
Dune::PDELab::Backend::Vector< GFS, typename P::range_type::field_type > range_type
The range type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:145
Overlapping parallel restarted GMRes solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:751
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
Y range_type
Definition: ovlpistlsolverbackend.hh:51
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: ovlpistlsolverbackend.hh:1031
Overlapping parallel BiCGStab solver with ILU0 preconditioner.
Definition: ovlpistlsolverbackend.hh:684
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:490
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlpistlsolverbackend.hh:898
virtual void apply(domain_type &v, const range_type &d)
Apply the preconditioner.
Definition: ovlpistlsolverbackend.hh:170
ISTLBackend_CG_AMG_SSOR(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Constructor.
Definition: ovlpistlsolverbackend.hh:1272
Definition: ovlpistlsolverbackend.hh:874
X domain_type
The domain type of the preconditioner.
Definition: ovlpistlsolverbackend.hh:205