dune-pdelab  2.7-git
assembledblockjacobipreconditioner.hh
Go to the documentation of this file.
1 // -*- tab-width: 2; indent-tabs-mode: nil -*-
2 // vi: set et ts=2 sw=2 sts=2:
3 #ifndef DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_MATRIXFREE_ASSEMBLEDBLOCKJACOBIPREONDITIONER_HH
5 
6 #if HAVE_EIGEN
7 
8 #include<dune/grid/common/scsgmapper.hh>
9 
11 
12 #include<Eigen/Dense>
13 
14 
15 namespace Dune {
16  namespace PDELab {
17 
18  namespace impl {
19  // Below we will assemble a block of the matrix as Eigen matrix. This
20  // accumulation view makes sure that we apply the weights.
21  template<typename C>
22  class WeightedEigenMatrixAccumulationView
23  {
24  public :
25  typedef C Container;
26 
27  typedef typename Container::Scalar value_type;
28  using size_type = typename Container::StorageIndex;
29  using weight_type = value_type;
30 
31  const weight_type& weight()
32  {
33  return _weight;
34  }
35 
36  auto data()
37  {
38  return _container.data();
39  }
40 
41  const auto data() const
42  {
43  return _container.data();
44  }
45 
46  auto& container()
47  {
48  return _container;
49  }
50 
51  template <typename LFSV, typename LFSU>
52  void accumulate(const LFSV& lfsv, size_type i, const LFSU& lfsu, size_type j, value_type value)
53  {
54  _container(i,j) += _weight * value;
55  }
56 
57  WeightedEigenMatrixAccumulationView(Container& container, weight_type weight)
58  : _container(container)
59  , _weight(weight)
60  {}
61 
62  private :
63  Container& _container;
64  weight_type _weight;
65  };
66  }
67 
68 
96  template<class BlockDiagonalLocalOperator, class GridFunctionSpace, class DomainField>
97  class AssembledBlockJacobiPreconditionerLocalOperator
100  {
101  // Extract some types
102  using GridView = typename GridFunctionSpace::Traits::GridViewType;
103  using Grid = typename GridView::Traits::Grid;
104  using EntityType = typename Grid::template Codim<0>::Entity;
105  using value_type = DomainField;
106  using DiagonalBlock = ::Eigen::Matrix<value_type, ::Eigen::Dynamic, ::Eigen::Dynamic, ::Eigen::RowMajor>;
107  using PermutationIndices = ::Eigen::Matrix<int, ::Eigen::Dynamic, 1>;
108  using TupleType = std::tuple<DiagonalBlock, PermutationIndices>;
109 
110  public :
111 
112  // Since this class implements something like D^{-1}v for a diagonal
113  // block matrix D we will only have volume methods. The underlying local
114  // operator that describes the block diagonal might of course have
115  // skeleton or boundary parts.
116  static constexpr bool doPatternVolume = true;
117  static constexpr bool doAlphaVolume = true;
118  static constexpr bool isLinear = BlockDiagonalLocalOperator::isLinear;
119 
130  AssembledBlockJacobiPreconditionerLocalOperator(const BlockDiagonalLocalOperator& blockDiagonalLocalOperator,
131  const GridFunctionSpace& gridFunctionSpace,
132  const bool verbose=0)
133  : _blockDiagonalLocalOperator(blockDiagonalLocalOperator)
134  , _gridFunctionSpace(gridFunctionSpace)
135  , _mapper(gridFunctionSpace.gridView())
136  , _precCache(_mapper.size())
137  , _verbose(verbose)
138  , _requireSetup(true)
139  {
140  }
141 
142 
143  // Before we can call the jacobian_apply methods we need to compute the
144  // LU decomposition of the diagonal block. We use the bool _requireSetup
145  // to make sure that we have computed this decomposition before we try to
146  // use it.
147  bool requireSetup()
148  {
149  return _requireSetup;
150  }
151  void setRequireSetup(bool requireSetup)
152  {
153  _requireSetup = requireSetup;
154  }
155 
156 
162  template<typename EG, typename LFSU, typename X, typename LFSV, typename Y>
163  void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, Y& y) const
164  {
165  // Matrix for storing the LU decomposition of the block diagonal
166  const int size = lfsu.size();
167  assert(lfsv.size() == size);
168 
169  // Assemble local block diagonal
170  DiagonalBlock diagonalBlock = DiagonalBlock::Constant(size, size, 0.0);
171  impl::WeightedEigenMatrixAccumulationView<DiagonalBlock> viewDiagonalBlock(diagonalBlock, y.weight());
172  assembleLocalDiagonalBlock(_blockDiagonalLocalOperator, eg, lfsu, x, lfsv, viewDiagonalBlock);
173 
174  // Compute the LU-factorization directly inside the matrix
175  ::Eigen::PartialPivLU<::Eigen::Ref<DiagonalBlock>> lu(diagonalBlock);
176  auto inversePermutationIndices = lu.permutationP().indices();
177 
178  // We need the inverse of the permutation matrix that Eigen gives
179  PermutationIndices permutationIndices(size);
180  for(int i=0; i<inversePermutationIndices.size(); ++i)
181  permutationIndices[inversePermutationIndices[i]] = i;
182 
183  // Fill up correct entry of preconditioner cache
184  std::size_t cache_idx = _mapper.index(eg.entity());
185  _precCache[cache_idx] = std::make_tuple(diagonalBlock, permutationIndices);
186  }
187 
188  template<typename EG, typename LFSU, typename X, typename Z, typename LFSV, typename Y>
189  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const X& x, const Z& z, const LFSV& lfsv, Y& y) const
190  {
191  // We don't need the linearization point as the preconditioner is
192  // already set up
193  jacobian_apply_volume(eg,lfsu,z,lfsv,y);
194  }
195 
196 
202  template<typename EG, typename LFSU, typename Z, typename LFSV, typename Y>
203  void jacobian_apply_volume(const EG& eg, const LFSU& lfsu, const Z& z, const LFSV& lfsv, Y& y) const
204  {
205  assert(not _requireSetup);
206 
207  // Get correct LU decomposition of the block diagonal from the cache
208  std::size_t cache_idx = _mapper.index(eg.entity());
209  const auto& cacheEntry = _precCache[cache_idx];
210  auto diagonalBlockLU = std::get<0>(cacheEntry);
211  auto permutationIndices = std::get<1>(cacheEntry);
212 
213  // Apply the preconditioner
214  const int size = lfsu.size();
215  const auto& z_container = accessBaseContainer(z);
216  std::vector<value_type> ztmp_container(size);
217  auto& y_container = accessBaseContainer(y);
218 
219  // Forward solve with lower triangle
220  for(int i=0; i<size; ++i){
221  value_type rhs(z_container[permutationIndices[i]]);
222  for(int j=0; j<i; ++j)
223  rhs -= diagonalBlockLU(i, j) * ztmp_container[j];
224  ztmp_container[i] = rhs; // L has ones on the diagonal
225  }
226  // Backward solve with upper triangular
227  for(int i=size-1; i>=0; i--){
228  value_type rhs(ztmp_container[i]);
229  for(int j=size-1; j>i; j--)
230  rhs -= diagonalBlockLU(i, j) * y_container[j];
231  y_container[i] = rhs/diagonalBlockLU(i,i);
232  }
233  }
234 
235  private :
236  BlockDiagonalLocalOperator _blockDiagonalLocalOperator;
237  const GridFunctionSpace& _gridFunctionSpace;
238  typename Dune::SingleCodimSingleGeomTypeMapper<GridView, 0> _mapper;
239  mutable std::vector<TupleType> _precCache;
240  const int _verbose;
241  bool _requireSetup;
242  };
243 
244  } // namespace PDELab
245 } // namespace Dune
246 
247 #endif // HAVE_EIGEN
248 
249 #endif
C & accessBaseContainer(C &c)
Definition: localvector.hh:368
GridView GridViewType
Definition: gridfunctionspace.hh:135
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void assembleLocalDiagonalBlock(const LocalOperator &localOperator, const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, MAT &mat)
A function for assembling a single diagonal block.
Definition: blockdiagonalwrapper.hh:226
Default flags for all local operators.
Definition: flags.hh:19
sparsity pattern generator
Definition: pattern.hh:14
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139