DUNE-FEM (unstable)

fastamgsmoother.hh
1// SPDX-FileCopyrightText: Copyright © 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
10namespace Dune
11{
12 namespace Amg
13 {
14
15 template<std::size_t level>
16 struct GaussSeidelPresmoothDefect {
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>
52 struct GaussSeidelPostsmoothDefect {
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
Dune namespace.
Definition: alignedallocator.hh:13
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)