dune-pdelab  2.7-git
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 #include <cstddef>
7 
8 #include <dune/common/deprecated.hh>
9 #include <dune/common/parallel/mpihelper.hh>
10 
11 #include <dune/grid/common/gridenums.hh>
12 
13 #include <dune/istl/io.hh>
14 #include <dune/istl/operators.hh>
15 #include <dune/istl/owneroverlapcopy.hh>
16 #include <dune/istl/paamg/amg.hh>
17 #include <dune/istl/paamg/pinfo.hh>
18 #include <dune/istl/preconditioners.hh>
19 #include <dune/istl/scalarproducts.hh>
20 #include <dune/istl/solvercategory.hh>
21 #include <dune/istl/solvers.hh>
22 #include <dune/istl/superlu.hh>
23 
31 
32 namespace Dune {
33  namespace PDELab {
34 
38 
39  //========================================================
40  // Generic support for nonoverlapping grids
41  //========================================================
42 
44 
52  template<typename GFS, typename M, typename X, typename Y>
54  : public Dune::AssembledLinearOperator<M,X,Y>
55  {
56  public:
64  typedef typename X::field_type field_type;
65 
67 
77  NonoverlappingOperator (const GFS& gfs_, const M& A)
78  : gfs(gfs_), _A_(A)
79  { }
80 
82 
87  virtual void apply (const X& x, Y& y) const override
88  {
89  using Backend::native;
90  // apply local operator; now we have sum y_p = sequential y
91  native(_A_).mv(native(x),native(y));
92 
93  // accumulate y on border
95  if (gfs.gridView().comm().size()>1)
96  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
97  }
98 
100 
105  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
106  {
107  using Backend::native;
108  // apply local operator; now we have sum y_p = sequential y
109  native(_A_).usmv(alpha,native(x),native(y));
110 
111  // accumulate y on border
113  if (gfs.gridView().comm().size()>1)
114  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
115  }
116 
117  SolverCategory::Category category() const override
118  {
119  return SolverCategory::nonoverlapping;
120  }
121 
123  virtual const M& getmat () const override
124  {
125  return _A_;
126  }
127 
128  private:
129  const GFS& gfs;
130  const M& _A_;
131  };
132 
133  // parallel scalar product assuming no overlap
134  template<class GFS, class X>
135  class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
136  {
137  public:
139  typedef X domain_type;
140  typedef typename X::ElementType field_type;
141 
142  SolverCategory::Category category() const override
143  {
144  return SolverCategory::nonoverlapping;
145  }
146 
149  NonoverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
150  : gfs(gfs_), helper(helper_)
151  {}
152 
157  virtual field_type dot (const X& x, const X& y) const override
158  {
159  // do local scalar product on unique partition
160  field_type sum = helper.disjointDot(x,y);
161 
162  // do global communication
163  return gfs.gridView().comm().sum(sum);
164  }
165 
169  virtual double norm (const X& x) const override
170  {
171  return sqrt(static_cast<double>(this->dot(x,x)));
172  }
173 
176  void make_consistent (X& x) const
177  {
179  if (gfs.gridView().comm().size()>1)
180  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
181  }
182 
183  private:
184  const GFS& gfs;
185  const ISTL::ParallelHelper<GFS>& helper;
186  };
187 
188  // parallel Richardson preconditioner
189  template<class GFS, class X, class Y>
190  class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
191  {
192  public:
194  typedef X domain_type;
196  typedef Y range_type;
198  typedef typename X::ElementType field_type;
199 
200  // define the category
201  SolverCategory::Category category() const override
202  {
203  return SolverCategory::nonoverlapping;
204  }
205 
207  NonoverlappingRichardson (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_)
208  : gfs(gfs_), helper(helper_)
209  {
210  }
211 
215  virtual void pre (X& x, Y& b) const override {}
216 
220  virtual void apply (X& v, const Y& d) const override
221  {
222  v = d;
223  }
224 
228  virtual void post (X& x) override {}
229 
230  private:
231  const GFS& gfs;
232  const ISTL::ParallelHelper<GFS>& helper;
233  };
234 
236 
249  template<typename A, typename X, typename Y>
251  : public Dune::Preconditioner<X,Y>
252  {
253 
255 
256  Diagonal _inverse_diagonal;
257 
258  public:
260 
264  typedef X domain_type;
266 
270  typedef Y range_type;
272  typedef typename X::ElementType field_type;
273 
274  SolverCategory::Category category() const override
275  {
276  return SolverCategory::nonoverlapping;
277  }
278 
280 
291  template<typename GFS>
292  NonoverlappingJacobi(const GFS& gfs, const A &m)
293  : _inverse_diagonal(m)
294  {
295  // make the diagonal consistent...
296  typename ISTL::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal);
297  gfs.gridView().communicate(addDH,
298  InteriorBorder_InteriorBorder_Interface,
299  ForwardCommunication);
300 
301  // ... and then invert it
302  _inverse_diagonal.invert();
303 
304  }
305 
307  virtual void pre (X& x, Y& b) override {}
308 
310  /*
311  * For this preconditioner, this method works with both consistent and
312  * inconsistent vectors: if d is consistent, v will be consistent, if d
313  * is inconsistent, v will be inconsistent.
314  */
315  virtual void apply (X& v, const Y& d) override
316  {
317  _inverse_diagonal.mv(d,v);
318  }
319 
321  virtual void post (X& x) override {}
322  };
323 
326 
328  template<class GFS>
330  {
332 
333  public:
340  explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
341  unsigned maxiter_=5000,
342  int verbose_=1)
343  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
344  {}
345 
350  template<class V>
351  typename V::ElementType norm (const V& v) const
352  {
353  V x(v); // make a copy because it has to be made consistent
355  PSP psp(gfs,phelper);
356  psp.make_consistent(x);
357  return psp.norm(x);
358  }
359 
367  template<class M, class V, class W>
368  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
369  {
371  POP pop(gfs,A);
373  PSP psp(gfs,phelper);
375  PRICH prich(gfs,phelper);
376  int verb=0;
377  if (gfs.gridView().comm().rank()==0) verb=verbose;
378  Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
379  Dune::InverseOperatorResult stat;
380  solver.apply(z,r,stat);
381  res.converged = stat.converged;
382  res.iterations = stat.iterations;
383  res.elapsed = stat.elapsed;
384  res.reduction = stat.reduction;
385  res.conv_rate = stat.conv_rate;
386  }
387 
390  {
391  return res;
392  }
393 
394  private:
395  const GFS& gfs;
396  PHELPER phelper;
398  unsigned maxiter;
399  int verbose;
400  };
401 
403  template<class GFS>
405  {
407 
408  const GFS& gfs;
409  PHELPER phelper;
411  unsigned maxiter;
412  int verbose;
413 
414  public:
416 
421  explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
422  unsigned maxiter_ = 5000,
423  int verbose_ = 1) :
424  gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
425  {}
426 
428 
433  template<class V>
434  typename V::ElementType norm (const V& v) const
435  {
436  V x(v); // make a copy because it has to be made consistent
438  PSP psp(gfs,phelper);
439  psp.make_consistent(x);
440  return psp.norm(x);
441  }
442 
444 
456  template<class M, class V, class W>
457  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
458  {
460  POP pop(gfs,A);
462  PSP psp(gfs,phelper);
463 
464  typedef NonoverlappingJacobi<M,V,W> PPre;
465  PPre ppre(gfs,Backend::native(A));
466 
467  int verb=0;
468  if (gfs.gridView().comm().rank()==0) verb=verbose;
469  CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
470  InverseOperatorResult stat;
471  solver.apply(z,r,stat);
472  res.converged = stat.converged;
473  res.iterations = stat.iterations;
474  res.elapsed = stat.elapsed;
475  res.reduction = stat.reduction;
476  res.conv_rate = stat.conv_rate;
477  }
478 
481  { return res; }
482  };
483 
485  template<class GFS>
487  {
489 
490  public:
497  explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
498  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
499  {}
500 
505  template<class V>
506  typename V::ElementType norm (const V& v) const
507  {
508  V x(v); // make a copy because it has to be made consistent
510  PSP psp(gfs,phelper);
511  psp.make_consistent(x);
512  return psp.norm(x);
513  }
514 
522  template<class M, class V, class W>
523  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
524  {
526  POP pop(gfs,A);
528  PSP psp(gfs,phelper);
530  PRICH prich(gfs,phelper);
531  int verb=0;
532  if (gfs.gridView().comm().rank()==0) verb=verbose;
533  Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
534  Dune::InverseOperatorResult stat;
535  solver.apply(z,r,stat);
536  res.converged = stat.converged;
537  res.iterations = stat.iterations;
538  res.elapsed = stat.elapsed;
539  res.reduction = stat.reduction;
540  res.conv_rate = stat.conv_rate;
541  }
542 
545  {
546  return res;
547  }
548 
549  private:
550  const GFS& gfs;
551  PHELPER phelper;
553  unsigned maxiter;
554  int verbose;
555  };
556 
557 
559  template<class GFS>
561  {
563 
564  public:
571  explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
572  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
573  {}
574 
579  template<class V>
580  typename V::ElementType norm (const V& v) const
581  {
582  V x(v); // make a copy because it has to be made consistent
584  PSP psp(gfs,phelper);
585  psp.make_consistent(x);
586  return psp.norm(x);
587  }
588 
596  template<class M, class V, class W>
597  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
598  {
600  POP pop(gfs,A);
602  PSP psp(gfs,phelper);
603 
604  typedef NonoverlappingJacobi<M,V,W> PPre;
605  PPre ppre(gfs,A);
606 
607  int verb=0;
608  if (gfs.gridView().comm().rank()==0) verb=verbose;
609  Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
610  Dune::InverseOperatorResult stat;
611  solver.apply(z,r,stat);
612  res.converged = stat.converged;
613  res.iterations = stat.iterations;
614  res.elapsed = stat.elapsed;
615  res.reduction = stat.reduction;
616  res.conv_rate = stat.conv_rate;
617  }
618 
621  {
622  return res;
623  }
624 
625  private:
626  const GFS& gfs;
627  PHELPER phelper;
629  unsigned maxiter;
630  int verbose;
631  };
632 
634  template<typename GFS>
636  {
638 
639  const GFS& gfs;
640  PHELPER phelper;
642 
643  public:
649  explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
650  : gfs(gfs_), phelper(gfs)
651  {}
652 
657  template<class V>
658  typename V::ElementType norm (const V& v) const
659  {
661  V x(v); // make a copy because it has to be made consistent
662  PSP psp(gfs,phelper);
663  psp.make_consistent(x);
664  return psp.norm(x);
665  }
666 
674  template<class M, class V, class W>
675  void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
676  {
677  Dune::SeqJac<M,V,W> jac(A,1,1.0);
678  jac.pre(z,r);
679  jac.apply(z,r);
680  jac.post(z);
681  if (gfs.gridView().comm().size()>1)
682  {
684  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
685  }
686  res.converged = true;
687  res.iterations = 1;
688  res.elapsed = 0.0;
689  res.reduction = reduction;
690  res.conv_rate = reduction; // pow(reduction,1.0/1)
691  }
692 
695  {
696  return res;
697  }
698  };
700 
701 
709  template<class GO,
710  template<class,class,class,int> class Preconditioner,
711  template<class> class Solver>
713  {
714  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
716 
717  public:
725  explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
726  : _grid_operator(grid_operator)
727  , gfs(grid_operator.trialGridFunctionSpace())
728  , phelper(gfs,verbose_)
729  , maxiter(maxiter_)
730  , steps(steps_)
731  , verbose(verbose_)
732  {}
733 
738  template<class Vector>
739  typename Vector::ElementType norm (const Vector& v) const
740  {
741  Vector x(v); // make a copy because it has to be made consistent
743  PSP psp(gfs,phelper);
744  psp.make_consistent(x);
745  return psp.norm(x);
746  }
747 
755  template<class M, class V, class W>
756  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
757  {
758  using MatrixType = Backend::Native<M>;
759  MatrixType& mat = Backend::native(A);
760  using VectorType = Backend::Native<W>;
761 #if HAVE_MPI
763  _grid_operator.make_consistent(A);
764  ISTL::assertParallelUG(gfs.gridView().comm());
765  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
766  phelper.createIndexSetAndProjectForAMG(mat, oocc);
767  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
768  Smoother smoother(mat, steps, 1.0);
769  typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP;
770  PSP psp(oocc);
771  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
772  Operator oop(mat,oocc);
773  typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother;
774  ParSmoother parsmoother(smoother, oocc);
775 #else
776  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
777  ParSmoother parsmoother(mat, steps, 1.0);
778  typedef Dune::SeqScalarProduct<VectorType> PSP;
779  PSP psp;
780  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
781  Operator oop(mat);
782 #endif
783  int verb=0;
784  if (gfs.gridView().comm().rank()==0) verb=verbose;
785  Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
786  Dune::InverseOperatorResult stat;
787  //make r consistent
788  if (gfs.gridView().comm().size()>1){
790  gfs.gridView().communicate(adddh,
791  Dune::InteriorBorder_InteriorBorder_Interface,
792  Dune::ForwardCommunication);
793  }
794 
795  solver.apply(z,r,stat);
796  res.converged = stat.converged;
797  res.iterations = stat.iterations;
798  res.elapsed = stat.elapsed;
799  res.reduction = stat.reduction;
800  res.conv_rate = stat.conv_rate;
801  }
802 
805  {
806  return res;
807  }
808 
809  private:
810  const GO& _grid_operator;
811  const GFS& gfs;
812  PHELPER phelper;
814  unsigned maxiter;
815  unsigned steps;
816  int verbose;
817  };
818 
821 
834  template<class GO>
836  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
837  {
838 
839  public:
847  explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
848  int steps_=5, int verbose_=1)
849  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_)
850  {}
851  };
852 
859  template<class GO>
861  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
862  {
863 
864  public:
872  explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
873  int steps_=5, int verbose_=1)
874  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_)
875  {}
876  };
879 
880  template<class GO,int s, template<class,class,class,int> class Preconditioner,
881  template<class> class Solver>
883  {
884  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
885  typedef typename ISTL::ParallelHelper<GFS> PHELPER;
886  typedef typename GO::Traits::Jacobian M;
887  typedef Backend::Native<M> MatrixType;
888  typedef typename GO::Traits::Domain V;
889  typedef Backend::Native<V> VectorType;
891 #if HAVE_MPI
892  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
893  typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother;
894  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
895 #else
896  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
897  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
898 #endif
899  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
900  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
901  typedef Dune::Amg::Parameters Parameters;
902 
903  public:
904  ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000,
905  int verbose_=1, bool reuse_=false,
906  bool usesuperlu_=true)
907  : _grid_operator(grid_operator)
908  , gfs(grid_operator.trialGridFunctionSpace())
909  , phelper(gfs,verbose_)
910  , maxiter(maxiter_)
911  , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu)
912  , verbose(verbose_)
913  , reuse(reuse_)
914  , firstapply(true)
915  , usesuperlu(usesuperlu_)
916  {
917  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
918  params.setDebugLevel(verbose_);
919 #if !HAVE_SUPERLU
920  if (phelper.rank() == 0 && usesuperlu == true)
921  {
922  std::cout << "WARNING: You are using AMG without SuperLU!"
923  << " Please consider installing SuperLU,"
924  << " or set the usesuperlu flag to false"
925  << " to suppress this warning." << std::endl;
926  }
927 #endif
928  }
929 
934  void setParameters(const Parameters& params_)
935  {
936  params = params_;
937  }
938 
946  const Parameters& parameters() const
947  {
948  return params;
949  }
950 
952  void setReuse(bool reuse_)
953  {
954  reuse = reuse_;
955  }
956 
958  bool getReuse() const
959  {
960  return reuse;
961  }
962 
963 
968  typename V::ElementType norm (const V& v) const
969  {
970  V x(v); // make a copy because it has to be made consistent
972  PSP psp(gfs,phelper);
973  psp.make_consistent(x);
974  return psp.norm(x);
975  }
976 
977  void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
978  {
979  Timer watch;
980  MatrixType& mat = Backend::native(A);
981  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
982  Dune::Amg::FirstDiagonal> > Criterion;
983 #if HAVE_MPI
984  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
985  _grid_operator.make_consistent(A);
986  phelper.createIndexSetAndProjectForAMG(A, oocc);
987  Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
988  Operator oop(mat, oocc);
989 #else
990  Comm oocc(gfs.gridView().comm());
991  Operator oop(mat);
992  Dune::SeqScalarProduct<VectorType> sp;
993 #endif
994  SmootherArgs smootherArgs;
995  smootherArgs.iterations = 1;
996  smootherArgs.relaxationFactor = 1;
997  //use noAccu or atOnceAccu
998  Criterion criterion(params);
999  stats.tprepare=watch.elapsed();
1000  watch.reset();
1001 
1002  int verb=0;
1003  if (gfs.gridView().comm().rank()==0) verb=verbose;
1004  //only construct a new AMG if the matrix changes
1005  if (reuse==false || firstapply==true){
1006  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
1007  firstapply = false;
1008  stats.tsetup = watch.elapsed();
1009  stats.levels = amg->maxlevels();
1010  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
1011  }
1012 
1013  Dune::InverseOperatorResult stat;
1014  // make r consistent
1015  if (gfs.gridView().comm().size()>1) {
1017  gfs.gridView().communicate(adddh,
1018  Dune::InteriorBorder_InteriorBorder_Interface,
1019  Dune::ForwardCommunication);
1020  }
1021  watch.reset();
1022  Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
1023  solver.apply(Backend::native(z),Backend::native(r),stat);
1024  stats.tsolve= watch.elapsed();
1025  res.converged = stat.converged;
1026  res.iterations = stat.iterations;
1027  res.elapsed = stat.elapsed;
1028  res.reduction = stat.reduction;
1029  res.conv_rate = stat.conv_rate;
1030  }
1031 
1037  {
1038  return stats;
1039  }
1040 
1041  private:
1042  const GO& _grid_operator;
1043  const GFS& gfs;
1044  PHELPER phelper;
1045  unsigned maxiter;
1046  Parameters params;
1047  int verbose;
1048  bool reuse;
1049  bool firstapply;
1050  bool usesuperlu;
1051  std::shared_ptr<AMG> amg;
1052  ISTLAMGStatistics stats;
1053  };
1054 
1057 
1072  template<class GO, int s=96>
1074  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1075  {
1076 
1077  public:
1078  ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1079  int verbose_=1, bool reuse_=false,
1080  bool usesuperlu_=true)
1081  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1082  {}
1083  };
1084 
1099  template<class GO, int s=96>
1101  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1102  {
1103 
1104  public:
1105  ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1106  int verbose_=1, bool reuse_=false,
1107  bool usesuperlu_=true)
1108  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1109  {}
1110  };
1111 
1126  template<class GO, int s=96>
1128  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
1129  {
1130 
1131  public:
1132  ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1133  int verbose_=1, bool reuse_=false,
1134  bool usesuperlu_=true)
1135  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1136  {}
1137  };
1140 
1141  } // namespace PDELab
1142 } // namespace Dune
1143 
1144 #endif // DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH
const std::string s
Definition: function.hh:843
void assertParallelUG(T comm)
Definition: parallelhelper.hh:445
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
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: blockmatrixdiagonal.hh:215
void invert()
Definition: blockmatrixdiagonal.hh:233
void mv(const X &x, Y &y) const
Definition: blockmatrixdiagonal.hh:239
Operator for the non-overlapping parallel case.
Definition: novlpistlsolverbackend.hh:55
NonoverlappingOperator(const GFS &gfs_, const M &A)
Construct a non-overlapping operator.
Definition: novlpistlsolverbackend.hh:77
virtual const M & getmat() const override
extract the matrix
Definition: novlpistlsolverbackend.hh:123
virtual void apply(const X &x, Y &y) const override
apply operator
Definition: novlpistlsolverbackend.hh:87
Backend::Native< X > domain_type
export type of vectors the matrix is applied to
Definition: novlpistlsolverbackend.hh:60
X::field_type field_type
export type of the entries for x
Definition: novlpistlsolverbackend.hh:64
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition: novlpistlsolverbackend.hh:105
Backend::Native< Y > range_type
export type of result vectors
Definition: novlpistlsolverbackend.hh:62
SolverCategory::Category category() const override
Definition: novlpistlsolverbackend.hh:117
Backend::Native< M > matrix_type
export type of matrix
Definition: novlpistlsolverbackend.hh:58
Definition: novlpistlsolverbackend.hh:136
X::ElementType field_type
Definition: novlpistlsolverbackend.hh:140
SolverCategory::Category category() const override
Definition: novlpistlsolverbackend.hh:142
X domain_type
export types
Definition: novlpistlsolverbackend.hh:139
NonoverlappingScalarProduct(const GFS &gfs_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: novlpistlsolverbackend.hh:149
virtual double norm(const X &x) const override
Norm of a right-hand side vector. The vector must be consistent on the interior+border partition.
Definition: novlpistlsolverbackend.hh:169
void make_consistent(X &x) const
make additive vector consistent
Definition: novlpistlsolverbackend.hh:176
virtual field_type dot(const X &x, const X &y) const override
Dot product of two vectors. It is assumed that the vectors are consistent on the interior+border part...
Definition: novlpistlsolverbackend.hh:157
Definition: novlpistlsolverbackend.hh:191
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:198
NonoverlappingRichardson(const GFS &gfs_, const ISTL::ParallelHelper< GFS > &helper_)
Constructor.
Definition: novlpistlsolverbackend.hh:207
virtual void apply(X &v, const Y &d) const override
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:220
X domain_type
The domain type of the preconditioner.
Definition: novlpistlsolverbackend.hh:194
Y range_type
The range type of the preconditioner.
Definition: novlpistlsolverbackend.hh:196
virtual void post(X &x) override
Clean up.
Definition: novlpistlsolverbackend.hh:228
virtual void pre(X &x, Y &b) const override
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:215
SolverCategory::Category category() const override
Definition: novlpistlsolverbackend.hh:201
parallel non-overlapping Jacobi preconditioner
Definition: novlpistlsolverbackend.hh:252
virtual void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:307
virtual void post(X &x) override
Clean up.
Definition: novlpistlsolverbackend.hh:321
virtual void apply(X &v, const Y &d) override
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:315
SolverCategory::Category category() const override
Definition: novlpistlsolverbackend.hh:274
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:272
X domain_type
The domain type of the operator.
Definition: novlpistlsolverbackend.hh:264
NonoverlappingJacobi(const GFS &gfs, const A &m)
Constructor.
Definition: novlpistlsolverbackend.hh:292
Y range_type
The range type of the operator.
Definition: novlpistlsolverbackend.hh:270
Nonoverlapping parallel CG solver without preconditioner.
Definition: novlpistlsolverbackend.hh:330
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:368
ISTLBackend_NOVLP_CG_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:340
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:389
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:351
Nonoverlapping parallel CG solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:405
const LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:480
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:434
ISTLBackend_NOVLP_CG_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:421
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:457
Nonoverlapping parallel BiCGStab solver without preconditioner.
Definition: novlpistlsolverbackend.hh:487
ISTLBackend_NOVLP_BCGS_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:497
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:523
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:544
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:506
Nonoverlapping parallel BiCGStab solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:561
ISTLBackend_NOVLP_BCGS_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:571
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:620
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:580
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:597
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: novlpistlsolverbackend.hh:636
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:694
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:658
ISTLBackend_NOVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: novlpistlsolverbackend.hh:649
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:675
Utility base class for preconditioned novlp backends.
Definition: novlpistlsolverbackend.hh:713
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:804
ISTLBackend_NOVLP_BASE_PREC(const GO &grid_operator, unsigned maxiter_=5000, unsigned steps_=5, int verbose_=1)
Constructor.
Definition: novlpistlsolverbackend.hh:725
Vector::ElementType norm(const Vector &v) const
Compute global norm of a vector.
Definition: novlpistlsolverbackend.hh:739
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
Solve the given linear system.
Definition: novlpistlsolverbackend.hh:756
Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:837
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:847
Nonoverlapping parallel CG solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:862
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:872
Definition: novlpistlsolverbackend.hh:883
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: novlpistlsolverbackend.hh:1036
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: novlpistlsolverbackend.hh:946
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
Definition: novlpistlsolverbackend.hh:977
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:968
void setParameters(const Parameters &params_)
set AMG parameters
Definition: novlpistlsolverbackend.hh:934
ISTLBackend_AMG_NOVLP(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:904
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: novlpistlsolverbackend.hh:958
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: novlpistlsolverbackend.hh:952
Nonoverlapping parallel CG solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1075
ISTLBackend_NOVLP_CG_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1078
Nonoverlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1102
ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1105
Nonoverlapping parallel LoopSolver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1129
ISTLBackend_NOVLP_LS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1132
Definition: parallelhelper.hh:51
RankIndex rank() const
Returns the MPI rank of this process.
Definition: parallelhelper.hh:248
void createIndexSetAndProjectForAMG(MatrixType &m, Comm &c)
Makes the matrix consistent and creates the parallel information for AMG.
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:429
Class providing some statistics of the AMG solver.
Definition: seqistlsolverbackend.hh:701
double tprepare
The needed for computing the parallel information and for adapting the linear system.
Definition: seqistlsolverbackend.hh:706
double tsetup
The time needed for building the AMG hierarchy (coarsening).
Definition: seqistlsolverbackend.hh:712
double tsolve
The time spent in solving the system (without building the hierarchy.
Definition: seqistlsolverbackend.hh:710
bool directCoarseLevelSolver
True if a direct solver was used on the coarset level.
Definition: seqistlsolverbackend.hh:716
int levels
the number of levels in the AMG hierarchy.
Definition: seqistlsolverbackend.hh:708
RFType conv_rate
Definition: solver.hh:36
bool converged
Definition: solver.hh:32
RFType reduction
Definition: solver.hh:35
unsigned int iterations
Definition: solver.hh:33
double elapsed
Definition: solver.hh:34
Definition: solver.hh:54
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:63
Definition: genericdatahandle.hh:667
Definition: recipe-operator-splitting.cc:108