dune-istl  3.0-git
fastamgsmoother.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_FASTAMGSMOOTHER_HH
4 #define DUNE_ISTL_FASTAMGSMOOTHER_HH
5 
6 #include <cstddef>
7 
8 namespace Dune
9 {
10  namespace Amg
11  {
12 
13  template<std::size_t level>
15 
16  template<typename M, typename X, typename Y>
17  static void apply(const M& A, X& x, Y& d,
18  const Y& b)
19  {
20  typedef typename M::ConstRowIterator RowIterator;
21  typedef typename M::ConstColIterator ColIterator;
22 
23  typename Y::iterator dIter=d.begin();
24  typename Y::const_iterator bIter=b.begin();
25  typename X::iterator xIter=x.begin();
26 
27  for(RowIterator row=A.begin(), end=A.end(); row != end;
28  ++row, ++dIter, ++xIter, ++bIter)
29  {
30  ColIterator col=(*row).begin();
31  *dIter = *bIter;
32 
33  for (; col.index()<row.index(); ++col)
34  (*col).mmv(x[col.index()],*dIter); // rhs -= sum_{j<i} a_ij * xnew_j
35  assert(row.index()==col.index());
36  ColIterator diag=col; // upper diagonal matrix not needed as x was 0 before.
37 
38  // Not recursive yet. Just solve with the diagonal
39  diag->solve(*xIter,*dIter);
40  *dIter=0; //as r=v
41 
42  // Update residual for the symmetric case
43  for(col=(*row).begin(); col.index()<row.index(); ++col)
44  col->mmv(*xIter, d[col.index()]); //d_j-=A_ij x_i
45  }
46  }
47  };
48 
49  template<std::size_t level>
51 
52  template<typename M, typename X, typename Y>
53  static void apply(const M& A, X& x, Y& d,
54  const Y& b)
55  {
56  typedef typename M::ConstRowIterator RowIterator;
57  typedef typename M::ConstColIterator ColIterator;
58  typedef typename Y::block_type YBlock;
59 
60  typename Y::iterator dIter=d.beforeEnd();
61  typename X::iterator xIter=x.beforeEnd();
62  typename Y::const_iterator bIter=b.beforeEnd();
63 
64  for(RowIterator row=A.beforeEnd(), end=A.beforeBegin(); row != end;
65  --row, --dIter, --xIter, --bIter)
66  {
67  ColIterator endCol=(*row).beforeBegin();
68  ColIterator col=(*row).beforeEnd();
69  *dIter = *bIter;
70 
71  for (; col.index()>row.index(); --col)
72  (*col).mmv(x[col.index()],*dIter); // rhs -= sum_{i>j} a_ij * xnew_j
73  assert(row.index()==col.index());
74  ColIterator diag=col;
75  YBlock v=*dIter;
76  // upper diagonal matrix
77  for (--col; col!=endCol; --col)
78  (*col).mmv(x[col.index()],v); // v -= sum_{j<i} a_ij * xold_j
79 
80  // Not recursive yet. Just solve with the diagonal
81  diag->solve(*xIter,v);
82 
83  *dIter-=v;
84 
85  // Update residual for the symmetric case
86  // Skip residual computation as it is not needed.
87  //for(col=(*row).begin();col.index()<row.index(); ++col)
88  //col.mmv(*xIter, d[col.index()]); //d_j-=A_ij x_i
89  }
90  }
91  };
92  } // end namespace Amg
93 } // end namespace Dune
94 #endif
Col col
Definition: matrixmatrix.hh:347
static void apply(const M &A, X &x, Y &d, const Y &b)
Definition: fastamgsmoother.hh:17
Definition: fastamgsmoother.hh:50
static void apply(const M &A, X &x, Y &d, const Y &b)
Definition: fastamgsmoother.hh:53
Definition: basearray.hh:19
Definition: fastamgsmoother.hh:14