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