dune-pdelab  2.0.0
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_NOVLPISTLSOLVERBACKEND_HH
4 #define DUNE_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 
32 
33 namespace Dune {
34  namespace PDELab {
35 
39 
40  //========================================================
41  // Generic support for nonoverlapping grids
42  //========================================================
43 
45 
53  template<typename GFS, typename M, typename X, typename Y>
55  : public Dune::AssembledLinearOperator<M,X,Y>
56  {
57  public:
59  typedef typename M::BaseT matrix_type;
61  typedef typename X::BaseT domain_type;
63  typedef typename Y::BaseT range_type;
65  typedef typename X::field_type field_type;
66 
67  //redefine the category, that is the only difference
68  enum {category=Dune::SolverCategory::nonoverlapping};
69 
71 
81  NonoverlappingOperator (const GFS& gfs_, const M& A)
82  : gfs(gfs_), _A_(A)
83  { }
84 
86 
91  virtual void apply (const X& x, Y& y) const
92  {
93  // apply local operator; now we have sum y_p = sequential y
94  istl::raw(_A_).mv(istl::raw(x),istl::raw(y));
95 
96  // accumulate y on border
98  if (gfs.gridView().comm().size()>1)
99  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
100  }
101 
103 
108  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
109  {
110  // apply local operator; now we have sum y_p = sequential y
111  istl::raw(_A_).usmv(alpha,istl::raw(x),istl::raw(y));
112 
113  // accumulate y on border
115  if (gfs.gridView().comm().size()>1)
116  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
117  }
118 
120  virtual const M& getmat () const
121  {
122  return _A_;
123  }
124 
125  private:
126  const GFS& gfs;
127  const M& _A_;
128  };
129 
130  // parallel scalar product assuming no overlap
131  template<class GFS, class X>
132  class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
133  {
134  public:
136  typedef X domain_type;
137  typedef typename X::ElementType field_type;
138 
140  enum {category=Dune::SolverCategory::nonoverlapping};
141 
144  NonoverlappingScalarProduct (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
145  : gfs(gfs_), helper(helper_)
146  {}
147 
152  virtual field_type dot (const X& x, const X& y)
153  {
154  // do local scalar product on unique partition
155  field_type sum = helper.disjointDot(x,y);
156 
157  // do global communication
158  return gfs.gridView().comm().sum(sum);
159  }
160 
164  virtual double norm (const X& x)
165  {
166  return sqrt(static_cast<double>(this->dot(x,x)));
167  }
168 
171  void make_consistent (X& x) const
172  {
174  if (gfs.gridView().comm().size()>1)
175  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
176  }
177 
178  private:
179  const GFS& gfs;
180  const istl::ParallelHelper<GFS>& helper;
181  };
182 
183  // parallel Richardson preconditioner
184  template<class GFS, class X, class Y>
185  class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
186  {
187  public:
189  typedef X domain_type;
191  typedef Y range_type;
193  typedef typename X::ElementType field_type;
194 
195  // define the category
196  enum {
198  category=Dune::SolverCategory::nonoverlapping
199  };
200 
202  NonoverlappingRichardson (const GFS& gfs_, const istl::ParallelHelper<GFS>& helper_)
203  : gfs(gfs_), helper(helper_)
204  {
205  }
206 
210  virtual void pre (X& x, Y& b) {}
211 
215  virtual void apply (X& v, const Y& d)
216  {
217  v = d;
218  }
219 
223  virtual void post (X& x) {}
224 
225  private:
226  const GFS& gfs;
227  const istl::ParallelHelper<GFS>& helper;
228  };
229 
231 
243  template<typename A, typename X, typename Y>
245  : public Dune::Preconditioner<X,Y>
246  {
247 
249 
250  Diagonal _inverse_diagonal;
251 
252  public:
254 
258  typedef X domain_type;
260 
264  typedef Y range_type;
266  typedef typename X::ElementType field_type;
267 
268  enum {
270  category=Dune::SolverCategory::nonoverlapping
271  };
272 
274 
285  template<typename GFS>
286  NonoverlappingJacobi(const GFS& gfs, const A &m)
287  : _inverse_diagonal(m)
288  {
289  // make the diagonal consistent...
290  typename istl::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal);
291  gfs.gridView().communicate(addDH,
292  InteriorBorder_InteriorBorder_Interface,
293  ForwardCommunication);
294 
295  // ... and then invert it
296  _inverse_diagonal.invert();
297 
298  }
299 
301  virtual void pre (X& x, Y& b) {}
302 
304  /*
305  * For this preconditioner, this method works with both consistent and
306  * inconsistent vectors: if d is consistent, v will be consistent, if d
307  * is inconsistent, v will be inconsistent.
308  */
309  virtual void apply (X& v, const Y& d)
310  {
311  _inverse_diagonal.mv(d,v);
312  }
313 
315  virtual void post (X& x) {}
316  };
317 
320 
322  template<class GFS>
324  {
326 
327  public:
334  explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
335  unsigned maxiter_=5000,
336  int verbose_=1)
337  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
338  {}
339 
344  template<class V>
345  typename V::ElementType norm (const V& v) const
346  {
347  V x(v); // make a copy because it has to be made consistent
349  PSP psp(gfs,phelper);
350  psp.make_consistent(x);
351  return psp.norm(x);
352  }
353 
361  template<class M, class V, class W>
362  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
363  {
365  POP pop(gfs,A);
367  PSP psp(gfs,phelper);
369  PRICH prich(gfs,phelper);
370  int verb=0;
371  if (gfs.gridView().comm().rank()==0) verb=verbose;
372  Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
373  Dune::InverseOperatorResult stat;
374  solver.apply(z,r,stat);
375  res.converged = stat.converged;
376  res.iterations = stat.iterations;
377  res.elapsed = stat.elapsed;
378  res.reduction = stat.reduction;
379  res.conv_rate = stat.conv_rate;
380  }
381 
384  {
385  return res;
386  }
387 
388  private:
389  const GFS& gfs;
390  PHELPER phelper;
392  unsigned maxiter;
393  int verbose;
394  };
395 
397  template<class GFS>
399  {
401 
402  const GFS& gfs;
403  PHELPER phelper;
405  unsigned maxiter;
406  int verbose;
407 
408  public:
410 
415  explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
416  unsigned maxiter_ = 5000,
417  int verbose_ = 1) :
418  gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
419  {}
420 
422 
427  template<class V>
428  typename V::ElementType norm (const V& v) const
429  {
430  V x(v); // make a copy because it has to be made consistent
432  PSP psp(gfs,phelper);
433  psp.make_consistent(x);
434  return psp.norm(x);
435  }
436 
438 
450  template<class M, class V, class W>
451  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
452  {
454  POP pop(gfs,A);
456  PSP psp(gfs,phelper);
457 
458  typedef NonoverlappingJacobi<M,V,W> PPre;
459  PPre ppre(gfs,istl::raw(A));
460 
461  int verb=0;
462  if (gfs.gridView().comm().rank()==0) verb=verbose;
463  CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
464  InverseOperatorResult stat;
465  solver.apply(z,r,stat);
466  res.converged = stat.converged;
467  res.iterations = stat.iterations;
468  res.elapsed = stat.elapsed;
469  res.reduction = stat.reduction;
470  res.conv_rate = stat.conv_rate;
471  }
472 
475  { return res; }
476  };
477 
479  template<class GFS>
481  {
483 
484  public:
491  explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
492  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
493  {}
494 
499  template<class V>
500  typename V::ElementType norm (const V& v) const
501  {
502  V x(v); // make a copy because it has to be made consistent
504  PSP psp(gfs,phelper);
505  psp.make_consistent(x);
506  return psp.norm(x);
507  }
508 
516  template<class M, class V, class W>
517  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
518  {
520  POP pop(gfs,A);
522  PSP psp(gfs,phelper);
524  PRICH prich(gfs,phelper);
525  int verb=0;
526  if (gfs.gridView().comm().rank()==0) verb=verbose;
527  Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
528  Dune::InverseOperatorResult stat;
529  solver.apply(z,r,stat);
530  res.converged = stat.converged;
531  res.iterations = stat.iterations;
532  res.elapsed = stat.elapsed;
533  res.reduction = stat.reduction;
534  res.conv_rate = stat.conv_rate;
535  }
536 
539  {
540  return res;
541  }
542 
543  private:
544  const GFS& gfs;
545  PHELPER phelper;
547  unsigned maxiter;
548  int verbose;
549  };
550 
551 
553  template<class GFS>
555  {
557 
558  public:
565  explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
566  : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
567  {}
568 
573  template<class V>
574  typename V::ElementType norm (const V& v) const
575  {
576  V x(v); // make a copy because it has to be made consistent
578  PSP psp(gfs,phelper);
579  psp.make_consistent(x);
580  return psp.norm(x);
581  }
582 
590  template<class M, class V, class W>
591  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
592  {
594  POP pop(gfs,A);
596  PSP psp(gfs,phelper);
597 
598  typedef NonoverlappingJacobi<M,V,W> PPre;
599  PPre ppre(gfs,A);
600 
601  int verb=0;
602  if (gfs.gridView().comm().rank()==0) verb=verbose;
603  Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
604  Dune::InverseOperatorResult stat;
605  solver.apply(z,r,stat);
606  res.converged = stat.converged;
607  res.iterations = stat.iterations;
608  res.elapsed = stat.elapsed;
609  res.reduction = stat.reduction;
610  res.conv_rate = stat.conv_rate;
611  }
612 
615  {
616  return res;
617  }
618 
619  private:
620  const GFS& gfs;
621  PHELPER phelper;
623  unsigned maxiter;
624  int verbose;
625  };
626 
628  template<typename GFS>
630  {
632 
633  const GFS& gfs;
634  PHELPER phelper;
636 
637  public:
643  explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
644  : gfs(gfs_), phelper(gfs)
645  {}
646 
651  template<class V>
652  typename V::ElementType norm (const V& v) const
653  {
655  V x(v); // make a copy because it has to be made consistent
656  PSP psp(gfs,phelper);
657  psp.make_consistent(x);
658  return psp.norm(x);
659  }
660 
668  template<class M, class V, class W>
669  void apply(M& A, V& z, W& r, typename W::ElementType reduction)
670  {
671  Dune::SeqJac<M,V,W> jac(A,1,1.0);
672  jac.pre(z,r);
673  jac.apply(z,r);
674  jac.post(z);
675  if (gfs.gridView().comm().size()>1)
676  {
678  gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
679  }
680  res.converged = true;
681  res.iterations = 1;
682  res.elapsed = 0.0;
683  res.reduction = reduction;
684  res.conv_rate = reduction; // pow(reduction,1.0/1)
685  }
686 
689  {
690  return res;
691  }
692  };
694 
695 
696  template<class GO,
697  template<class,class,class,int> class Preconditioner,
698  template<class> class Solver>
700  {
701  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
703 
704  public:
712  explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
713  : _grid_operator(grid_operator)
714  , gfs(grid_operator.trialGridFunctionSpace())
715  , phelper(gfs,verbose_)
716  , maxiter(maxiter_)
717  , steps(steps_)
718  , verbose(verbose_)
719  {}
720 
725  template<class Vector>
726  typename Vector::ElementType norm (const Vector& v) const
727  {
728  Vector x(v); // make a copy because it has to be made consistent
730  PSP psp(gfs,phelper);
731  psp.make_consistent(x);
732  return psp.norm(x);
733  }
734 
742  template<class M, class V, class W>
743  void apply(M& A, V& z, W& r, typename V::ElementType reduction)
744  {
745  typedef typename M::BaseT MatrixType;
746  MatrixType& mat=istl::raw(A);
747  typedef typename V::BaseT VectorType;
748 #if HAVE_MPI
750  _grid_operator.make_consistent(A);
751  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
752  phelper.createIndexSetAndProjectForAMG(mat, oocc);
753  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
754  Smoother smoother(mat, steps, 1.0);
755  typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP;
756  PSP psp(oocc);
757  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
758  Operator oop(mat,oocc);
759  typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother;
760  ParSmoother parsmoother(smoother, oocc);
761 #else
762  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
763  ParSmoother parsmoother(mat, steps, 1.0);
764  typedef Dune::SeqScalarProduct<VectorType> PSP;
765  PSP psp;
766  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
767  Operator oop(mat);
768 #endif
769  int verb=0;
770  if (gfs.gridView().comm().rank()==0) verb=verbose;
771  Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
772  Dune::InverseOperatorResult stat;
773  //make r consistent
774  if (gfs.gridView().comm().size()>1){
776  gfs.gridView().communicate(adddh,
777  Dune::InteriorBorder_InteriorBorder_Interface,
778  Dune::ForwardCommunication);
779  }
780 
781  solver.apply(z,r,stat);
782  res.converged = stat.converged;
783  res.iterations = stat.iterations;
784  res.elapsed = stat.elapsed;
785  res.reduction = stat.reduction;
786  res.conv_rate = stat.conv_rate;
787  }
788 
791  {
792  return res;
793  }
794 
795  private:
796  const GO& _grid_operator;
797  const GFS& gfs;
798  PHELPER phelper;
800  unsigned maxiter;
801  unsigned steps;
802  int verbose;
803  };
804 
807 
819  template<class GO>
821  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
822  {
823 
824  public:
832  explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
833  int steps_=5, int verbose_=1)
834  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_)
835  {}
836  };
837 
841  template<class GO>
843  : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
844  {
845 
846  public:
854  explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000,
855  int steps_=5, int verbose_=1)
856  : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_)
857  {}
858  };
861 
862  template<class GO,int s, template<class,class,class,int> class Preconditioner,
863  template<class> class Solver>
865  {
866  typedef typename GO::Traits::TrialGridFunctionSpace GFS;
867  typedef typename istl::ParallelHelper<GFS> PHELPER;
868  typedef typename GO::Traits::Jacobian M;
869  typedef typename M::BaseT MatrixType;
870  typedef typename GO::Traits::Domain V;
871  typedef typename V::BaseT VectorType;
873 #if HAVE_MPI
874  typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
875  typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother;
876  typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
877 #else
878  typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
879  typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
880 #endif
881  typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
882  typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
883  typedef Dune::Amg::Parameters Parameters;
884 
885  public:
886  ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000,
887  int verbose_=1, bool reuse_=false,
888  bool usesuperlu_=true)
889  : _grid_operator(grid_operator)
890  , gfs(grid_operator.trialGridFunctionSpace())
891  , phelper(gfs,verbose_)
892  , maxiter(maxiter_)
893  , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu)
894  , verbose(verbose_)
895  , reuse(reuse_)
896  , firstapply(true)
897  , usesuperlu(usesuperlu_)
898  {
899  params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
900  params.setDebugLevel(verbose_);
901 #if !HAVE_SUPERLU
902  if (phelper.rank() == 0 && usesuperlu == true)
903  {
904  std::cout << "WARNING: You are using AMG without SuperLU!"
905  << " Please consider installing SuperLU,"
906  << " or set the usesuperlu flag to false"
907  << " to suppress this warning." << std::endl;
908  }
909 #endif
910  }
911 
916  void setParameters(const Parameters& params_)
917  {
918  params = params_;
919  }
920 
921  void setparams(Parameters params_) DUNE_DEPRECATED_MSG("setparams() is deprecated, use setParameters() instead")
922  {
923  params = params_;
924  }
925 
933  const Parameters& parameters() const
934  {
935  return params;
936  }
937 
942  typename V::ElementType norm (const V& v) const
943  {
944  V x(v); // make a copy because it has to be made consistent
946  PSP psp(gfs,phelper);
947  psp.make_consistent(x);
948  return psp.norm(x);
949  }
950 
951  void apply(M& A, V& z, V& r, typename V::ElementType reduction)
952  {
953  Timer watch;
954  MatrixType& mat=istl::raw(A);
955  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
956  Dune::Amg::FirstDiagonal> > Criterion;
957 #if HAVE_MPI
958  Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
959  _grid_operator.make_consistent(A);
960  phelper.createIndexSetAndProjectForAMG(A, oocc);
961  Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
962  Operator oop(mat, oocc);
963 #else
964  Comm oocc(gfs.gridView().comm());
965  Operator oop(mat);
966  Dune::SeqScalarProduct<VectorType> sp;
967 #endif
968  SmootherArgs smootherArgs;
969  smootherArgs.iterations = 1;
970  smootherArgs.relaxationFactor = 1;
971  //use noAccu or atOnceAccu
972  Criterion criterion(params);
973  stats.tprepare=watch.elapsed();
974  watch.reset();
975 
976  int verb=0;
977  if (gfs.gridView().comm().rank()==0) verb=verbose;
978  //only construct a new AMG if the matrix changes
979  if (reuse==false || firstapply==true){
980  amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
981  firstapply = false;
982  stats.tsetup = watch.elapsed();
983  stats.levels = amg->maxlevels();
984  stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
985  }
986 
987  Dune::InverseOperatorResult stat;
988  // make r consistent
989  if (gfs.gridView().comm().size()>1) {
991  gfs.gridView().communicate(adddh,
992  Dune::InteriorBorder_InteriorBorder_Interface,
993  Dune::ForwardCommunication);
994  }
995  watch.reset();
996  Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
997  solver.apply(istl::raw(z),istl::raw(r),stat);
998  stats.tsolve= watch.elapsed();
999  res.converged = stat.converged;
1000  res.iterations = stat.iterations;
1001  res.elapsed = stat.elapsed;
1002  res.reduction = stat.reduction;
1003  res.conv_rate = stat.conv_rate;
1004  }
1005 
1011  {
1012  return stats;
1013  }
1014 
1015  private:
1016  const GO& _grid_operator;
1017  const GFS& gfs;
1018  PHELPER phelper;
1019  unsigned maxiter;
1020  Parameters params;
1021  int verbose;
1022  bool reuse;
1023  bool firstapply;
1024  bool usesuperlu;
1025  Dune::shared_ptr<AMG> amg;
1026  ISTLAMGStatistics stats;
1027  };
1028 
1031 
1044  template<class GO, int s=96>
1046  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
1047  {
1048 
1049  public:
1050  ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1051  int verbose_=1, bool reuse_=false,
1052  bool usesuperlu_=true)
1053  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1054  {}
1055  };
1056 
1069  template<class GO, int s=96>
1071  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
1072  {
1073 
1074  public:
1075  ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1076  int verbose_=1, bool reuse_=false,
1077  bool usesuperlu_=true)
1078  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1079  {}
1080  };
1081 
1094  template<class GO, int s=96>
1096  : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
1097  {
1098 
1099  public:
1100  ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000,
1101  int verbose_=1, bool reuse_=false,
1102  bool usesuperlu_=true)
1103  : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_)
1104  {}
1105  };
1108 
1109  } // namespace PDELab
1110 } // namespace Dune
1111 
1112 #endif
NonoverlappingRichardson(const GFS &gfs_, const istl::ParallelHelper< GFS > &helper_)
Constructor.
Definition: novlpistlsolverbackend.hh:202
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
Solve the given linear system.
Definition: novlpistlsolverbackend.hh:743
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:614
Nonoverlapping parallel CG solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:398
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:451
X domain_type
The domain type of the operator.
Definition: novlpistlsolverbackend.hh:258
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const
apply operator to x, scale and add:
Definition: novlpistlsolverbackend.hh:108
Class providing some statistics of the AMG solver.
Definition: seqistlsolverbackend.hh:472
void mv(const X &x, Y &y) const
Definition: blockmatrixdiagonal.hh:202
M::BaseT matrix_type
export type of matrix
Definition: novlpistlsolverbackend.hh:59
double tsetup
The time needed for building the AMG hierarchy (coarsening).
Definition: seqistlsolverbackend.hh:484
Solver to be used for explicit time-steppers with (block-)diagonal mass matrix.
Definition: novlpistlsolverbackend.hh:629
void make_consistent(X &x) const
make additive vector consistent
Definition: novlpistlsolverbackend.hh:171
bool converged
Definition: solver.hh:32
virtual void apply(const X &x, Y &y) const
apply operator
Definition: novlpistlsolverbackend.hh:91
const LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:474
double tsolve
The time spent in solving the system (without building the hierarchy.
Definition: seqistlsolverbackend.hh:482
ISTLBackend_NOVLP_LS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1100
virtual void post(X &x)
Clean up.
Definition: novlpistlsolverbackend.hh:223
NonoverlappingScalarProduct(const GFS &gfs_, const istl::ParallelHelper< GFS > &helper_)
Constructor needs to know the grid function space.
Definition: novlpistlsolverbackend.hh:144
X domain_type
export types
Definition: novlpistlsolverbackend.hh:136
ISTLBackend_NOVLP_BASE_PREC(const GO &grid_operator, unsigned maxiter_=5000, unsigned steps_=5, int verbose_=1)
Constructor.
Definition: novlpistlsolverbackend.hh:712
Nonoverlapping parallel BiCGStab solver with Jacobi preconditioner.
Definition: novlpistlsolverbackend.hh:554
double elapsed
Definition: solver.hh:34
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:538
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:215
int levels
the number of levels in the AMG hierarchy.
Definition: seqistlsolverbackend.hh:480
Definition: novlpistlsolverbackend.hh:132
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:152
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:164
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:832
X::field_type field_type
export type of the entries for x
Definition: novlpistlsolverbackend.hh:65
virtual const M & getmat() const
extract the matrix
Definition: novlpistlsolverbackend.hh:120
void apply(M &A, V &z, W &r, typename W::ElementType reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:669
NonoverlappingJacobi(const GFS &gfs, const A &m)
Constructor.
Definition: novlpistlsolverbackend.hh:286
Nonoverlapping parallel LoopSolver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1095
ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1075
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:210
const ISTLAMGStatistics & statistics() const
Get statistics of the AMG solver (no of levels, timings).
Definition: novlpistlsolverbackend.hh:1010
NonoverlappingOperator(const GFS &gfs_, const M &A)
Construct a non-overlapping operator.
Definition: novlpistlsolverbackend.hh:81
Nonoverlapping parallel CG solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:842
ISTLBackend_NOVLP_BCGS_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:491
void setparams(Parameters params_)
Definition: novlpistlsolverbackend.hh:921
Definition: novlpistlsolverbackend.hh:68
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:942
The category the preconditioner is part of.
Definition: novlpistlsolverbackend.hh:198
Nonoverlapping parallel CG solver preconditioned with AMG smoothed by SSOR.
Definition: novlpistlsolverbackend.hh:1045
V & raw(V &v)
Returns the raw ISTL object associated with v, or v itself it is already an ISTL object.
Definition: backend/istl/utility.hh:26
ISTLBackend_NOVLP_ExplicitDiagonal(const GFS &gfs_)
make a linear solver object
Definition: novlpistlsolverbackend.hh:643
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52
void setParameters(const Parameters &params_)
set AMG parameters
Definition: novlpistlsolverbackend.hh:916
X domain_type
The domain type of the preconditioner.
Definition: novlpistlsolverbackend.hh:189
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:266
Y range_type
The range type of the operator.
Definition: novlpistlsolverbackend.hh:264
ISTLBackend_NOVLP_CG_NOPREC(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:334
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:854
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:383
Definition: novlpistlsolverbackend.hh:864
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:591
void invert()
Definition: blockmatrixdiagonal.hh:196
parallel non-overlapping Jacobi preconditioner
Definition: novlpistlsolverbackend.hh:244
Nonoverlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR. ...
Definition: novlpistlsolverbackend.hh:1070
Definition: parallelhelper.hh:45
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:428
Definition: genericdatahandle.hh:623
virtual void post(X &x)
Clean up.
Definition: novlpistlsolverbackend.hh:315
void createIndexSetAndProjectForAMG(MatrixType &m, Comm &c)
Makes the matrix consistent and creates the parallel information for AMG.
Y::BaseT range_type
export type of result vectors
Definition: novlpistlsolverbackend.hh:63
RFType reduction
Definition: solver.hh:35
Definition: novlpistlsolverbackend.hh:185
Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR.
Definition: novlpistlsolverbackend.hh:820
X::ElementType field_type
Definition: novlpistlsolverbackend.hh:137
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:652
ISTLBackend_NOVLP_CG_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:415
ISTLBackend_NOVLP_BCGS_Jacobi(const GFS &gfs_, unsigned maxiter_=5000, int verbose_=1)
make a linear solver object
Definition: novlpistlsolverbackend.hh:565
Definition: solver.hh:42
Nonoverlapping parallel CG solver without preconditioner.
Definition: novlpistlsolverbackend.hh:323
The category the preconditioner is part of.
Definition: novlpistlsolverbackend.hh:270
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:500
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:790
virtual void apply(X &v, const Y &d)
Apply the precondioner.
Definition: novlpistlsolverbackend.hh:309
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:345
bool directCoarseLevelSolver
True if a direct solver was used on the coarset level.
Definition: seqistlsolverbackend.hh:488
X::ElementType field_type
The field type of the preconditioner.
Definition: novlpistlsolverbackend.hh:193
Definition: novlpistlsolverbackend.hh:699
Vector::ElementType norm(const Vector &v) const
Compute global norm of a vector.
Definition: novlpistlsolverbackend.hh:726
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:517
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: novlpistlsolverbackend.hh:301
const Dune::PDELab::LinearSolverResult< double > & result() const
Return access to result data.
Definition: novlpistlsolverbackend.hh:688
Definition: novlpistlsolverbackend.hh:140
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
Definition: blockmatrixdiagonal.hh:177
RankIndex rank() const
Returns the MPI rank of this process.
Definition: parallelhelper.hh:226
double tprepare
The needed for computing the parallel information and for adapting the linear system.
Definition: seqistlsolverbackend.hh:478
ISTLBackend_NOVLP_CG_AMG_SSOR(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:1050
Y range_type
The range type of the preconditioner.
Definition: novlpistlsolverbackend.hh:191
X::BaseT domain_type
export type of vectors the matrix is applied to
Definition: novlpistlsolverbackend.hh:61
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: novlpistlsolverbackend.hh:933
RFType conv_rate
Definition: solver.hh:36
ISTLBackend_AMG_NOVLP(const GO &grid_operator, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: novlpistlsolverbackend.hh:886
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
Definition: novlpistlsolverbackend.hh:951
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: novlpistlsolverbackend.hh:574
void apply(M &A, V &z, W &r, typename V::ElementType reduction)
solve the given linear system
Definition: novlpistlsolverbackend.hh:362
Operator for the non-overlapping parallel case.
Definition: novlpistlsolverbackend.hh:54
const std::string s
Definition: function.hh:1103
unsigned int iterations
Definition: solver.hh:33
Nonoverlapping parallel BiCGStab solver without preconditioner.
Definition: novlpistlsolverbackend.hh:480