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