dune-pdelab  2.5-dev
novlpistlsolverbackend.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_NOVLPISTLSOLVERBACKEND_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH
5 
6 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
7 #include "ensureistlinclude.hh"
8 
9 #include <cstddef>
10 
11 #include <dune/common/deprecated.hh>
12 #include <dune/common/parallel/mpihelper.hh>
13 
14 #include <dune/grid/common/gridenums.hh>
15 
16 #include <dune/istl/io.hh>
17 #include <dune/istl/operators.hh>
18 #include <dune/istl/owneroverlapcopy.hh>
19 #include <dune/istl/paamg/amg.hh>
20 #include <dune/istl/paamg/pinfo.hh>
21 #include <dune/istl/preconditioners.hh>
22 #include <dune/istl/scalarproducts.hh>
23 #include <dune/istl/solvercategory.hh>
24 #include <dune/istl/solvers.hh>
25 #include <dune/istl/superlu.hh>
26 
34 
35 namespace Dune {
36  namespace PDELab {
37 
41 
42  //========================================================
43  // Generic support for nonoverlapping grids
44  //========================================================
45 
47 
55  template<typename GFS, typename M, typename X, typename Y>
57  : public Dune::AssembledLinearOperator<M,X,Y>
58  {
59  public:
67  typedef typename X::field_type field_type;
68 
69  //redefine the category, that is the only difference
70  enum {category=Dune::SolverCategory::nonoverlapping};
71 
73 
83  NonoverlappingOperator (const GFS& gfs_, const M& A)
84  : gfs(gfs_), _A_(A)
85  { }
86 
88 
93  virtual void apply (const X& x, Y& y) const
94  {
95  using Backend::native;
96  // apply local operator; now we have sum y_p = sequential y
97  native(_A_).mv(native(x),native(y));
98 
99  // accumulate y on border
101  if (gfs.gridView().comm().size()>1)
102  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
103  }
104 
106 
111  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
112  {
113  using Backend::native;
114  // apply local operator; now we have sum y_p = sequential y
115  native(_A_).usmv(alpha,native(x),native(y));
116 
117  // accumulate y on border
119  if (gfs.gridView().comm().size()>1)
120  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
121  }
122 
124  virtual const M& getmat () const
125  {
126  return _A_;
127  }
128 
129  private:
130  const GFS& gfs;
131  const M& _A_;
132  };
133 
134  // parallel scalar product assuming no overlap
135  template<class GFS, class X>
136  class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
137  {
138  public:
140  typedef X domain_type;
141  typedef typename X::ElementType field_type;
142 
144  enum {category=Dune::SolverCategory::nonoverlapping};
145 
148  NonoverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
149  : gfs(gfs_), helper(helper_)
150  {}
151 
156  virtual field_type dot (const X& x, const X& y)
157  {
158  // do local scalar product on unique partition
159  field_type sum = helper.disjointDot(x,y);
160 
161  // do global communication
162  return gfs.gridView().comm().sum(sum);
163  }
164 
168  virtual double norm (const X& x)
169  {
170  return sqrt(static_cast<double>(this->dot(x,x)));
171  }
172 
175  void make_consistent (X& x) const
176  {
178  if (gfs.gridView().comm().size()>1)
179  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
180  }
181 
182  private:
183  const GFS& gfs;
184  const ISTL::ParallelHelper<GFS>& helper;
185  };
186 
187  // parallel Richardson preconditioner
188  template<class GFS, class X, class Y>
189  class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
190  {
191  public:
193  typedef X domain_type;
195  typedef Y range_type;
197  typedef typename X::ElementType field_type;
198 
199  // define the category
200  enum {
202  category=Dune::SolverCategory::nonoverlapping
203  };
204 
206  NonoverlappingRichardson (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
207  : gfs(gfs_), helper(helper_)
208  {
209  }
210 
214  virtual void pre (X& x, Y& b) {}
215 
219  virtual void apply (X& v, const Y& d)
220  {
221  v = d;
222  }
223 
227  virtual void post (X& x) {}
228 
229  private:
230  const GFS& gfs;
231  const ISTL::ParallelHelper<GFS>& helper;
232  };
233 
235 
248  template<typename A, typename X, typename Y>
250  : public Dune::Preconditioner<X,Y>
251  {
252 
254 
255  Diagonal _inverse_diagonal;
256 
257  public:
259 
263  typedef X domain_type;
265 
269  typedef Y range_type;
271  typedef typename X::ElementType field_type;
272 
273  enum {
275  category=Dune::SolverCategory::nonoverlapping
276  };
277 
279 
290  template<typename GFS>
291  NonoverlappingJacobi(const GFS& gfs, const A &m)
292  : _inverse_diagonal(m)
293  {
294  // make the diagonal consistent...
295  typename ISTL::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal);
296  gfs.gridView().communicate(addDH,
297  InteriorBorder_InteriorBorder_Interface,
298  ForwardCommunication);
299 
300  // ... and then invert it
301  _inverse_diagonal.invert();
302 
303  }
304 
306  virtual void pre (X& x, Y& b) {}
307 
309  /*
310  * For this preconditioner, this method works with both consistent and
311  * inconsistent vectors: if d is consistent, v will be consistent, if d
312  * is inconsistent, v will be inconsistent.
313  */
314  virtual void apply (X& v, const Y& d)
315  {
316  _inverse_diagonal.mv(d,v);
317  }
318 
320  virtual void post (X& x) {}
321  };
322 
325 
327  template<class GFS>
329  {
331 
332  public:
339  explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
340  unsigned maxiter_=5000,
341  int verbose_=1)
342  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
343  {}
344 
349  template<class V>
350  typename V::ElementType norm (const V& v) const
351  {
352  V x(v); // make a copy because it has to be made consistent
354  PSP psp(gfs,phelper);
355  psp.make_consistent(x);
356  return psp.norm(x);
357  }
358 
366  template<class M, class V, class W>
367  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
368  {
370  POP pop(gfs,A);
372  PSP psp(gfs,phelper);
374  PRICH prich(gfs,phelper);
375  int verb=0;
376  if (gfs.gridView().comm().rank()==0) verb=verbose;
377  Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
378  Dune::InverseOperatorResult stat;
379  solver.apply(z,r,stat);
380  res.converged = stat.converged;
381  res.iterations = stat.iterations;
382  res.elapsed = stat.elapsed;
383  res.reduction = stat.reduction;
384  res.conv_rate = stat.conv_rate;
385  }
386 
389  {
390  return res;
391  }
392 
393  private:
394  const GFS& gfs;
395  PHELPER phelper;
397  unsigned maxiter;
398  int verbose;
399  };
400 
402  template<class GFS>
404  {
406 
407  const GFS& gfs;
408  PHELPER phelper;
410  unsigned maxiter;
411  int verbose;
412 
413  public:
415 
420  explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
421  unsigned maxiter_ = 5000,
422  int verbose_ = 1) :
423  gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
424  {}
425 
427 
432  template<class V>
433  typename V::ElementType norm (const V& v) const
434  {
435  V x(v); // make a copy because it has to be made consistent
437  PSP psp(gfs,phelper);
438  psp.make_consistent(x);
439  return psp.norm(x);
440  }
441 
443 
455  template<class M, class V, class W>
456  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
457  {
459  POP pop(gfs,A);
461  PSP psp(gfs,phelper);
462 
463  typedef NonoverlappingJacobi<M,V,W> PPre;
464  PPre ppre(gfs,Backend::native(A));
465 
466  int verb=0;
467  if (gfs.gridView().comm().rank()==0) verb=verbose;
468  CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
469  InverseOperatorResult stat;
470  solver.apply(z,r,stat);
471  res.converged = stat.converged;
472  res.iterations = stat.iterations;
473  res.elapsed = stat.elapsed;
474  res.reduction = stat.reduction;
475  res.conv_rate = stat.conv_rate;
476  }
477 
480  { return res; }
481  };
482 
484  template<class GFS>
486  {
488 
489  public:
496  explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
497  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
498  {}
499 
504  template<class V>
505  typename V::ElementType norm (const V& v) const
506  {
507  V x(v); // make a copy because it has to be made consistent
509  PSP psp(gfs,phelper);
510  psp.make_consistent(x);
511  return psp.norm(x);
512  }
513 
521  template<class M, class V, class W>
522  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
523  {
525  POP pop(gfs,A);
527  PSP psp(gfs,phelper);
529  PRICH prich(gfs,phelper);
530  int verb=0;
531  if (gfs.gridView().comm().rank()==0) verb=verbose;
532  Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
533  Dune::InverseOperatorResult stat;
534  solver.apply(z,r,stat);
535  res.converged = stat.converged;
536  res.iterations = stat.iterations;
537  res.elapsed = stat.elapsed;
538  res.reduction = stat.reduction;
539  res.conv_rate = stat.conv_rate;
540  }
541 
544  {
545  return res;
546  }
547 
548  private:
549  const GFS& gfs;
550  PHELPER phelper;
552  unsigned maxiter;
553  int verbose;
554  };
555 
556 
558  template<class GFS>
560  {
562 
563  public:
570  explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
571  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
572  {}
573 
578  template<class V>
579  typename V::ElementType norm (const V& v) const
580  {
581  V x(v); // make a copy because it has to be made consistent
583  PSP psp(gfs,phelper);
584  psp.make_consistent(x);
585  return psp.norm(x);
586  }
587 
595  template<class M, class V, class W>
596  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
597  {
599  POP pop(gfs,A);
601  PSP psp(gfs,phelper);
602 
603  typedef NonoverlappingJacobi<M,V,W> PPre;
604  PPre ppre(gfs,A);
605 
606  int verb=0;
607  if (gfs.gridView().comm().rank()==0) verb=verbose;
608  Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
609  Dune::InverseOperatorResult stat;
610  solver.apply(z,r,stat);
611  res.converged = stat.converged;
612  res.iterations = stat.iterations;
613  res.elapsed = stat.elapsed;
614  res.reduction = stat.reduction;
615  res.conv_rate = stat.conv_rate;
616  }
617 
620  {
621  return res;
622  }
623 
624  private:
625  const GFS& gfs;
626  PHELPER phelper;
628  unsigned maxiter;
629  int verbose;
630  };
631 
633  template<typename GFS>
635  {
637 
638  const GFS& gfs;
639  PHELPER phelper;
641 
642  public:
648  explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
649  : gfs(gfs_), phelper(gfs)
650  {}
651 
656  template<class V>
657  typename V::ElementType norm (const V& v) const
658  {
660  V x(v); // make a copy because it has to be made consistent
661  PSP psp(gfs,phelper);
662  psp.make_consistent(x);
663  return psp.norm(x);
664  }
665 
673  template<class M, class V, class W>
674  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
675  {
676  Dune::SeqJac<M,V,W> jac(A,1,1.0);
677  jac.pre(z,r);
678  jac.apply(z,r);
679  jac.post(z);
680  if (gfs.gridView().comm().size()>1)
681  {
683  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
684  }
685  res.converged = true;
686  res.iterations = 1;
687  res.elapsed = 0.0;
688  res.reduction = reduction;
689  res.conv_rate = reduction; // pow(reduction,1.0/1)
690  }
691 
694  {
695  return res;
696  }
697  };
699 
700 
708  template<class GO,
709  template<class,class,class,int> class Preconditioner,
710  template<class> class Solver>
712  {
713  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
715 
716  public:
724  explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
725  : _grid_operator(grid_operator)
726  , gfs(grid_operator.trialGridFunctionSpace())
727  , phelper(gfs,verbose_)
728  , maxiter(maxiter_)
729  , steps(steps_)
730  , verbose(verbose_)
731  {}
732 
737  template<class Vector>
738  typename Vector::ElementType norm (const Vector& v) const
739  {
740  Vector x(v); // make a copy because it has to be made consistent
742  PSP psp(gfs,phelper);
743  psp.make_consistent(x);
744  return psp.norm(x);
745  }
746 
754  template<class M, class V, class W>
755  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
756  {
757  using MatrixType = Backend::Native<M>;
758  MatrixType& mat = Backend::native(A);
759  using VectorType = Backend::Native<W>;
760 #if HAVE_MPI
762  _grid_operator.make_consistent(A);
763  ISTL::assertParallelUG(gfs.gridView().comm());
764  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
765  phelper.createIndexSetAndProjectForAMG(mat, oocc);
766  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
767  Smoother smoother(mat, steps, 1.0);
768  typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP;
769  PSP psp(oocc);
770  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
771  Operator oop(mat,oocc);
772  typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother;
773  ParSmoother parsmoother(smoother, oocc);
774 #else
775  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
776  ParSmoother parsmoother(mat, steps, 1.0);
777  typedef Dune::SeqScalarProduct<VectorType> PSP;
778  PSP psp;
779  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
780  Operator oop(mat);
781 #endif
782  int verb=0;
783  if (gfs.gridView().comm().rank()==0) verb=verbose;
784  Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
785  Dune::InverseOperatorResult stat;
786  //make r consistent
787  if (gfs.gridView().comm().size()>1){
789  gfs.gridView().communicate(adddh,
790  Dune::InteriorBorder_InteriorBorder_Interface,
791  Dune::ForwardCommunication);
792  }
793 
794  solver.apply(z,r,stat);
795  res.converged = stat.converged;
796  res.iterations = stat.iterations;
797  res.elapsed = stat.elapsed;
798  res.reduction = stat.reduction;
799  res.conv_rate = stat.conv_rate;
800  }
801 
804  {
805  return res;
806  }
807 
808  private:
809  const GO& _grid_operator;
810  const GFS& gfs;
811  PHELPER phelper;
813  unsigned maxiter;
814  unsigned steps;
815  int verbose;
816  };
817 
820 
833  template<class GO>
835  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
836  {
837 
838  public:
846  explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
847  int steps_=5, int verbose_=1)
848  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_)
849  {}
850  };
851 
858  template<class GO>
860  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
861  {
862 
863  public:
871  explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
872  int steps_=5, int verbose_=1)
873  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_)
874  {}
875  };
878 
879  template<class GO,int s, template<class,class,class,int> class Preconditioner,
880  template<class> class Solver>
882  {
883  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
884  typedef typename ISTL::ParallelHelper<GFS> PHELPER;
885  typedef typename GO::Traits::Jacobian M;
886  typedef Backend::Native<M> MatrixType;
887  typedef typename GO::Traits::Domain V;
888  typedef Backend::Native<V> VectorType;
890 #if HAVE_MPI
891  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
892  typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother;
893  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
894 #else
895  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
896  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
897 #endif
898  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
899  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
900  typedef Dune::Amg::Parameters Parameters;
901 
902  public:
903  ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000,
904  int verbose_=1, bool reuse_=false,
905  bool usesuperlu_=true)
906  : _grid_operator(grid_operator)
907  , gfs(grid_operator.trialGridFunctionSpace())
908  , phelper(gfs,verbose_)
909  , maxiter(maxiter_)
910  , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu)
911  , verbose(verbose_)
912  , reuse(reuse_)
913  , firstapply(true)
914  , usesuperlu(usesuperlu_)
915  {
916  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
917  params.setDebugLevel(verbose_);
918 #if !HAVE_SUPERLU
919  if (phelper.rank() == 0 && usesuperlu == true)
920  {
921  std::cout << "WARNING: You are using AMG without SuperLU!"
922  << " Please consider installing SuperLU,"
923  << " or set the usesuperlu flag to false"
924  << " to suppress this warning." << std::endl;
925  }
926 #endif
927  }
928 
933  void setParameters(const Parameters& params_)
934  {
935  params = params_;
936  }
937 
945  const Parameters& parameters() const
946  {
947  return params;
948  }
949 
951  void setReuse(bool reuse_)
952  {
953  reuse = reuse_;
954  }
955 
957  bool getReuse() const
958  {
959  return reuse;
960  }
961 
962 
967  typename V::ElementType norm (const V& v) const
968  {
969  V x(v); // make a copy because it has to be made consistent
971  PSP psp(gfs,phelper);
972  psp.make_consistent(x);
973  return psp.norm(x);
974  }
975 
976  void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
977  {
978  Timer watch;
979  MatrixType& mat = Backend::native(A);
980  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
981  Dune::Amg::FirstDiagonal> > Criterion;
982 #if HAVE_MPI
983  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
984  _grid_operator.make_consistent(A);
985  phelper.createIndexSetAndProjectForAMG(A, oocc);
986  Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
987  Operator oop(mat, oocc);
988 #else
989  Comm oocc(gfs.gridView().comm());
990  Operator oop(mat);
991  Dune::SeqScalarProduct<VectorType> sp;
992 #endif
993  SmootherArgs smootherArgs;
994  smootherArgs.iterations = 1;
995  smootherArgs.relaxationFactor = 1;
996  //use noAccu or atOnceAccu
997  Criterion criterion(params);
998  stats.tprepare=watch.elapsed();
999  watch.reset();
1000 
1001  int verb=0;
1002  if (gfs.gridView().comm().rank()==0) verb=verbose;
1003  //only construct a new AMG if the matrix changes
1004  if (reuse==false || firstapply==true){
1005  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1006  firstapply = false;
1007  stats.tsetup = watch.elapsed();
1008  stats.levels = amg->maxlevels();
1009  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1010  }
1011 
1012  Dune::InverseOperatorResult stat;
1013  // make r consistent
1014  if (gfs.gridView().comm().size()>1) {
1016  gfs.gridView().communicate(adddh,
1017  Dune::InteriorBorder_InteriorBorder_Interface,
1018  Dune::ForwardCommunication);
1019  }
1020  watch.reset();
1021  Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
1022  solver.apply(Backend::native(z),Backend::native(r),stat);
1023  stats.tsolve= watch.elapsed();
1024  res.converged = stat.converged;
1025  res.iterations = stat.iterations;
1026  res.elapsed = stat.elapsed;
1027  res.reduction = stat.reduction;
1028  res.conv_rate = stat.conv_rate;
1029  }
1030 
1036  {
1037  return stats;
1038  }
1039 
1040  private:
1041  const GO& _grid_operator;
1042  const GFS& gfs;
1043  PHELPER phelper;
1044  unsigned maxiter;
1045  Parameters params;
1046  int verbose;
1047  bool reuse;
1048  bool firstapply;
1049  bool usesuperlu;
1050  std::shared_ptr<AMG> amg;
1051  ISTLAMGStatistics stats;
1052  };
1053 
1056 
1071  template<class GO, int s=96>
1073  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1074  {
1075 
1076  public:
1077  ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1078  int verbose_=1, bool reuse_=false,
1079  bool usesuperlu_=true)
1080  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1081  {}
1082  };
1083 
1098  template<class GO, int s=96>
1100  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1101  {
1102 
1103  public:
1104  ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1105  int verbose_=1, bool reuse_=false,
1106  bool usesuperlu_=true)
1107  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1108  {}
1109  };
1110 
1125  template<class GO, int s=96>
1127  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
1128  {
1129 
1130  public:
1131  ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1132  int verbose_=1, bool reuse_=false,
1133  bool usesuperlu_=true)
1134  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1135  {}
1136  };
1139 
1140  } // namespace PDELab
1141 } // namespace Dune
1142 
1143 #endif // DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:803
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:388
parallel non-overlapping Jacobi preconditioner
Definition: novlpistlsolverbackend.hh:249
Operator for the non-overlapping parallel case.
Definition: novlpistlsolverbackend.hh:56
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:433
NonoverlappingRichardson(const GFS &gfs_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor.
Definition: novlpistlsolverbackend.hh:206
ISTLBackend_NOVLP_BASE_PREC(const GO &grid_operator, unsigned maxiter_=5000, unsigned steps_=5, int verbose_=1)
Constructor.
Definition: novlpistlsolverbackend.hh:724
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:214
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:314
double elapsed
Definition: solver.hh:34
const std::string s
Definition: function.hh:830
X::field_type field_type
export type of the entries for x
Definition: novlpistlsolverbackend.hh:67
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:306
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:219
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
Definition: novlpistlsolverbackend.hh:189
ISTLBackend_NOVLP_CG_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:420
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:197
Backend::Native< M > matrix_type
export type of matrix
Definition: novlpistlsolverbackend.hh:61
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:693
void assertParallelUG(T comm)
Definition: parallelhelper.hh:433
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpistlsolverbackend.hh:111
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
unsigned int iterations
Definition: solver.hh:33
Nonoverlapping parallel CG solver without preconditioner.
Definition: novlpistlsolverbackend.hh:328
X domain_type
The domain type of the operator.
Definition: novlpistlsolverbackend.hh:263
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: novlpistlsolverbackend.hh:945
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: novlpistlsolverbackend.hh:957
ISTLBackend_NOVLP_CG_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:339
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename W::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:674
ISTLBackend_NOVLP_CG_SSORk(const GO &grid_operator, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:871
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
Definition: novlpistlsolverbackend.hh:976
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
Solve the given linear system.
Definition: novlpistlsolverbackend.hh:755
Y range_type
The range type of the operator.
Definition: novlpistlsolverbackend.hh:269
Y range_type
The range type of the preconditioner.
Definition: novlpistlsolverbackend.hh:195
ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1104
Utility base class for preconditioned novlp backends.
Definition: novlpistlsolverbackend.hh:711
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: novlpistlsolverbackend.hh:951
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:350
NonoverlappingOperator(const GFS &gfs_, const M &A)
Construct a non-overlapping operator.
Definition: novlpistlsolverbackend.hh:83
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:619
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:505
NonoverlappingJacobi(const GFS &gfs, const A &m)
Constructor.
Definition: novlpistlsolverbackend.hh:291
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:967
Definition: blockmatrixdiagonal.hh:217
const LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:479
Nonoverlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR. ...
Definition: novlpistlsolverbackend.hh:1099
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:456
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:271
Backend::Native< Y > range_type
export type of result vectors
Definition: novlpistlsolverbackend.hh:65
Backend::Native< X > domain_type
export type of vectors the matrix is applied to
Definition: novlpistlsolverbackend.hh:63
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: novlpistlsolverbackend.hh:156
Nonoverlapping parallel BiCGStab solver without preconditioner.
Definition: novlpistlsolverbackend.hh:485
virtual void post(X &x)
Clean up.
Definition: novlpistlsolverbackend.hh:320
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:367
RFType conv_rate
Definition: solver.hh:36
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: novlpistlsolverbackend.hh:1035
Vector::ElementType norm(const Vector &v) const
Compute global norm of a vector.
Definition: novlpistlsolverbackend.hh:738
void setParameters(const Parameters &params_)
set AMG parameters
Definition: novlpistlsolverbackend.hh:933
Nonoverlapping parallel CG solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:859
ISTLBackend_NOVLP_BCGS_SSORk(const GO &grid_operator, unsigned maxiter_=5000, int steps_=5, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:846
virtual void apply(const X &x, Y &y) const
apply operator
Definition: novlpistlsolverbackend.hh:93
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: novlpistlsolverbackend.hh:634
ISTLBackend_NOVLP_CG_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1077
Nonoverlapping parallel CG solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:403
ISTLBackend_NOVLP_BCGS_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:496
void invert()
Definition: blockmatrixdiagonal.hh:236
virtual const M & getmat() const
extract the matrix
Definition: novlpistlsolverbackend.hh:124
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:522
virtual void post(X &x)
Clean up.
Definition: novlpistlsolverbackend.hh:227
Nonoverlapping parallel CG solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1072
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:657
Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:834
ISTLBackend_NOVLP_BCGS_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:570
virtual double norm(const X &x)
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition...
Definition: novlpistlsolverbackend.hh:168
Nonoverlapping parallel BiCGStab solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:559
Definition: novlpistlsolverbackend.hh:70
bool converged
Definition: solver.hh:32
X domain_type
The domain type of the preconditioner.
Definition: novlpistlsolverbackend.hh:193
X::ElementType field_type
Definition: novlpistlsolverbackend.hh:141
ISTLBackend_NOVLP_LS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1131
ISTLBackend_NOVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: novlpistlsolverbackend.hh:648
void make_consistent(X &x) const
make additive vector consistent
Definition: novlpistlsolverbackend.hh:175
Class providing some statistics of the AMG solver.
Definition: seqistlsolverbackend.hh:544
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:543
Nonoverlapping parallel LoopSolver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1126
void apply(M &A, V &z, W &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:596
Definition: genericdatahandle.hh:622
Definition: novlpistlsolverbackend.hh:881
void mv(const X &x, Y &y) const
Definition: blockmatrixdiagonal.hh:242
Definition: novlpistlsolverbackend.hh:136
RFType reduction
Definition: solver.hh:35
Definition: solver.hh:42
X domain_type
export types
Definition: novlpistlsolverbackend.hh:140
ISTLBackend_AMG_NOVLP(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:903
Definition: parallelhelper.hh:53
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:579
NonoverlappingScalarProduct(const GFS &gfs_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: novlpistlsolverbackend.hh:148