dune-istl  2.2.1
preconditioners.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_PRECONDITIONERS_HH
4 #define DUNE_PRECONDITIONERS_HH
5 
6 #include<cmath>
7 #include<complex>
8 #include<iostream>
9 #include<iomanip>
10 #include<string>
11 
12 #include"solvercategory.hh"
13 #include "istlexception.hh"
14 #include "matrixutils.hh"
15 #include "io.hh"
16 #include "gsetc.hh"
17 #include "ilu.hh"
18 
19 
20 namespace Dune {
58  //=====================================================================
72  //=====================================================================
73  template<class X, class Y>
75  public:
77  typedef X domain_type;
79  typedef Y range_type;
81  typedef typename X::field_type field_type;
82 
98  virtual void pre (X& x, Y& b) = 0;
99 
110  virtual void apply (X& v, const Y& d) = 0;
111 
120  virtual void post (X& x) = 0;
121 
122  // every abstract base class has a virtual destructor
123  virtual ~Preconditioner () {}
124  };
125 
126 
127  //=====================================================================
128  // Implementation of this interface for sequential ISTL-preconditioners
129  //=====================================================================
130 
131 
143  template<class M, class X, class Y, int l=1>
144  class SeqSSOR : public Preconditioner<X,Y> {
145  public:
147  typedef M matrix_type;
149  typedef X domain_type;
151  typedef Y range_type;
153  typedef typename X::field_type field_type;
154 
155  // define the category
156  enum {
159 
167  SeqSSOR (const M& A, int n, field_type w)
168  : _A_(A), _n(n), _w(w)
169  {
171  }
172 
178  virtual void pre (X& x, Y& b) {}
179 
185  virtual void apply (X& v, const Y& d)
186  {
187  for (int i=0; i<_n; i++){
188  bsorf(_A_,v,d,_w,BL<l>());
189  bsorb(_A_,v,d,_w,BL<l>());
190  }
191  }
192 
198  virtual void post (X& x) {}
199 
200  private:
202  const M& _A_;
204  int _n;
206  field_type _w;
207  };
208 
209 
210 
222  template<class M, class X, class Y, int l=1>
223  class SeqSOR : public Preconditioner<X,Y> {
224  public:
226  typedef M matrix_type;
228  typedef X domain_type;
230  typedef Y range_type;
232  typedef typename X::field_type field_type;
233 
234  // define the category
235  enum {
238  };
239 
247  SeqSOR (const M& A, int n, field_type w)
248  : _A_(A), _n(n), _w(w)
249  {
251  }
252 
258  virtual void pre (X& x, Y& b) {}
259 
265  virtual void apply (X& v, const Y& d)
266  {
267  this->template apply<true>(v,d);
268  }
269 
278  template<bool forward>
279  void apply(X& v, const Y& d)
280  {
281  if(forward)
282  for (int i=0; i<_n; i++){
283  bsorf(_A_,v,d,_w,BL<l>());
284  }
285  else
286  for (int i=0; i<_n; i++){
287  bsorb(_A_,v,d,_w,BL<l>());
288  }
289  }
290 
296  virtual void post (X& x) {}
297 
298  private:
300  const M& _A_;
302  int _n;
304  field_type _w;
305  };
306 
307 
318  template<class M, class X, class Y, int l=1>
319  class SeqGS : public Preconditioner<X,Y> {
320  public:
322  typedef M matrix_type;
324  typedef X domain_type;
326  typedef Y range_type;
328  typedef typename X::field_type field_type;
329 
330  // define the category
331  enum {
334  };
335 
343  SeqGS (const M& A, int n, field_type w)
344  : _A_(A), _n(n), _w(w)
345  {
347  }
348 
354  virtual void pre (X& x, Y& b) {}
355 
361  virtual void apply (X& v, const Y& d)
362  {
363  for (int i=0; i<_n; i++){
364  dbgs(_A_,v,d,_w,BL<l>());
365  }
366  }
367 
373  virtual void post (X& x) {}
374 
375  private:
377  const M& _A_;
379  int _n;
381  field_type _w;
382  };
383 
384 
395  template<class M, class X, class Y, int l=1>
396  class SeqJac : public Preconditioner<X,Y> {
397  public:
399  typedef M matrix_type;
401  typedef X domain_type;
403  typedef Y range_type;
405  typedef typename X::field_type field_type;
406 
407  // define the category
408  enum {
411  };
412 
420  SeqJac (const M& A, int n, field_type w)
421  : _A_(A), _n(n), _w(w)
422  {
424  }
425 
431  virtual void pre (X& x, Y& b) {}
432 
438  virtual void apply (X& v, const Y& d)
439  {
440  for (int i=0; i<_n; i++){
441  dbjac(_A_,v,d,_w,BL<l>());
442  }
443  }
444 
450  virtual void post (X& x) {}
451 
452  private:
454  const M& _A_;
456  int _n;
458  field_type _w;
459  };
460 
461 
462 
474  template<class M, class X, class Y, int l=1>
475  class SeqILU0 : public Preconditioner<X,Y> {
476  public:
480  typedef X domain_type;
482  typedef Y range_type;
484  typedef typename X::field_type field_type;
485 
486  // define the category
487  enum {
490  };
491 
498  SeqILU0 (const M& A, field_type w)
499  : ILU(A) // copy A
500  {
501  _w =w;
502  bilu0_decomposition(ILU);
503  }
504 
510  virtual void pre (X& x, Y& b) {}
511 
517  virtual void apply (X& v, const Y& d)
518  {
519  bilu_backsolve(ILU,v,d);
520  v *= _w;
521  }
522 
528  virtual void post (X& x) {}
529 
530  private:
532  field_type _w;
534  matrix_type ILU;
535  };
536 
537 
551  template<class M, class X, class Y, int l=1>
552  class SeqILUn : public Preconditioner<X,Y> {
553  public:
557  typedef X domain_type;
559  typedef Y range_type;
561  typedef typename X::field_type field_type;
562 
563  // define the category
564  enum {
567  };
568 
576  SeqILUn (const M& A, int n, field_type w)
577  : ILU(A.N(),A.M(),M::row_wise)
578  {
579  _n = n;
580  _w = w;
581  bilu_decomposition(A,n,ILU);
582  }
583 
589  virtual void pre (X& x, Y& b) {}
590 
596  virtual void apply (X& v, const Y& d)
597  {
598  bilu_backsolve(ILU,v,d);
599  v *= _w;
600  }
601 
607  virtual void post (X& x) {}
608 
609  private:
611  matrix_type ILU;
613  int _n;
615  field_type _w;
616  };
617 
618 
619 
628  template<class X, class Y>
630  public:
632  typedef X domain_type;
634  typedef Y range_type;
636  typedef typename X::field_type field_type;
637 
638  // define the category
639  enum {
642  };
643 
649  Richardson (field_type w=1.0)
650  {
651  _w = w;
652  }
653 
659  virtual void pre (X& x, Y& b) {}
660 
666  virtual void apply (X& v, const Y& d)
667  {
668  v = d;
669  v *= _w;
670  }
671 
677  virtual void post (X& x) {}
678 
679  private:
681  field_type _w;
682  };
683 
684 
685 
688 } // end namespace
689 
690 #endif