Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
schur_matrix.h
1 // ---------------------------------------------------------------------
2 // @f$Id: schur_matrix.h 30748 2013-09-16 21:34:46Z bangerth @f$
3 //
4 // Copyright (C) 2001 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__schur_matrix_h
18 #define __deal2__schur_matrix_h
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/subscriptor.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/base/logstream.h>
24 #include <deal.II/lac/vector_memory.h>
25 #include <deal.II/lac/block_vector.h>
26 #include <vector>
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 
97 template <class MA_inverse, class MB, class MDt, class MC>
98 class SchurMatrix : public Subscriptor
99 {
100 public:
101 
115  SchurMatrix(const MA_inverse &Ainv,
116  const MB &B,
117  const MDt &Dt,
118  const MC &C,
120  const std::vector<unsigned int> &signature = std::vector<unsigned int>(0));
121 
140  void prepare_rhs (BlockVector<double> &dst,
141  const BlockVector<double> &src) const;
142 
147  void vmult (BlockVector<double> &dst,
148  const BlockVector<double> &src) const;
149 
150 // void Tmult(BlockVector<double>& dst, const BlockVector<double>& src) const;
151 
156  double residual (BlockVector<double> &dst,
157  const BlockVector<double> &src,
158  const BlockVector<double> &rhs) const;
159 
165  void postprocess (BlockVector<double> &dst,
166  const BlockVector<double> &src,
167  const BlockVector<double> &rhs) const;
168 
177  void debug_level(unsigned int l);
178 private:
187 
208 
212  std::vector<types::global_dof_index> signature;
213 
217  unsigned int debug;
218 };
219 
221 //---------------------------------------------------------------------------
222 
223 template <class MA_inverse, class MB, class MDt, class MC>
225 ::SchurMatrix(const MA_inverse &Ainv,
226  const MB &B,
227  const MDt &Dt,
228  const MC &C,
230  const std::vector<unsigned int> &signature)
231  : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C),
232  mem(mem),
233  signature(signature),
234  debug(0)
235 {
236 }
237 
238 
239 template <class MA_inverse, class MB, class MDt, class MC>
240 void
242 ::debug_level(unsigned int l)
243 {
244  debug = l;
245 }
246 
247 
248 template <class MA_inverse, class MB, class MDt, class MC>
251  const BlockVector<double> &src) const
252 {
253  deallog.push("Schur");
254  if (debug > 0)
255  deallog << "src:" << src.l2_norm() << std::endl;
256 
257  C->vmult(dst, src);
258  if (debug > 0)
259  deallog << "C:" << dst.l2_norm() << std::endl;
260 
261  BlockVector<double> *h1 = mem.alloc();
262  if (signature.size()>0)
263  h1->reinit(signature);
264  else
265  h1->reinit(B->n_block_cols(), src.block(0).size());
266  Dt->Tvmult(*h1,src);
267  if (debug > 0)
268  deallog << "Dt:" << h1->l2_norm() << std::endl;
269 
270  BlockVector<double> *h2 = mem.alloc();
271  h2->reinit(*h1);
272  Ainv->vmult(*h2, *h1);
273  if (debug > 0)
274  deallog << "Ainverse:" << h2->l2_norm() << std::endl;
275 
276  mem.free(h1);
277  B->vmult_add(dst, *h2);
278  if (debug > 0)
279  deallog << "dst:" << dst.l2_norm() << std::endl;
280 
281  mem.free(h2);
282  deallog.pop();
283 }
284 
285 
286 template <class MA_inverse, class MB, class MDt, class MC>
289  const BlockVector<double> &src,
290  const BlockVector<double> &rhs) const
291 {
292  vmult(dst, src);
293  dst.scale(-1.);
294  dst += rhs;
295  return dst.l2_norm();
296 }
297 
298 
299 template <class MA_inverse, class MB, class MDt, class MC>
302  const BlockVector<double> &src) const
303 {
304  Assert (src.n_blocks() == B->n_block_cols(),
305  ExcDimensionMismatch(src.n_blocks(), B->n_block_cols()));
306  Assert (dst.n_blocks() == B->n_block_rows(),
307  ExcDimensionMismatch(dst.n_blocks(), B->n_block_rows()));
308 
309  deallog.push("Schur-prepare");
310  if (debug > 0)
311  deallog << "src:" << src.l2_norm() << std::endl;
312  BlockVector<double> *h1 = mem.alloc();
313  if (signature.size()>0)
314  h1->reinit(signature);
315  else
316  h1->reinit(B->n_block_cols(), src.block(0).size());
317  Ainv->vmult(*h1, src);
318  if (debug > 0)
319  deallog << "Ainverse:" << h1->l2_norm() << std::endl;
320  B->vmult_add(dst, *h1);
321  if (debug > 0)
322  deallog << "dst:" << dst.l2_norm() << std::endl;
323  mem.free(h1);
324  deallog.pop();
325 }
326 
327 
328 template <class MA_inverse, class MB, class MDt, class MC>
331  const BlockVector<double> &src,
332  const BlockVector<double> &rhs) const
333 {
334  Assert (dst.n_blocks() == B->n_block_cols(),
335  ExcDimensionMismatch(dst.n_blocks(), B->n_block_cols()));
336  Assert (rhs.n_blocks() == B->n_block_cols(),
337  ExcDimensionMismatch(rhs.n_blocks(), B->n_block_cols()));
338  Assert (src.n_blocks() == B->n_block_rows(),
339  ExcDimensionMismatch(src.n_blocks(), B->n_block_rows()));
340 
341  deallog.push("Schur-post");
342  if (debug > 0)
343  deallog << "src:" << src.l2_norm() << std::endl;
344  BlockVector<double> *h1 = mem.alloc();
345  if (signature.size()>0)
346  h1->reinit(signature);
347  else
348  h1->reinit(B->n_block_cols(), src.block(0).size());
349  Dt->Tvmult(*h1, src);
350  if (debug > 0)
351  deallog << "Dt:" << h1->l2_norm() << std::endl;
352  h1->sadd(-1.,rhs);
353  Ainv->vmult(dst,*h1);
354  if (debug > 0)
355  deallog << "dst:" << dst.l2_norm() << std::endl;
356  mem.free(h1);
357  deallog.pop();
358 }
359 
360 
361 DEAL_II_NAMESPACE_CLOSE
362 
363 #endif
void pop()
void vmult(BlockVector< double > &dst, const BlockVector< double > &src) const
Definition: schur_matrix.h:250
std::vector< types::global_dof_index > signature
Definition: schur_matrix.h:212
void postprocess(BlockVector< double > &dst, const BlockVector< double > &src, const BlockVector< double > &rhs) const
Definition: schur_matrix.h:330
void prepare_rhs(BlockVector< double > &dst, const BlockVector< double > &src) const
Definition: schur_matrix.h:301
double residual(BlockVector< double > &dst, const BlockVector< double > &src, const BlockVector< double > &rhs) const
Definition: schur_matrix.h:288
SchurMatrix(const MA_inverse &Ainv, const MB &B, const MDt &Dt, const MC &C, VectorMemory< BlockVector< double > > &mem, const std::vector< unsigned int > &signature=std::vector< unsigned int >(0))
Definition: schur_matrix.h:225
unsigned int debug
Definition: schur_matrix.h:217
const SmartPointer< const MB, SchurMatrix< MA_inverse, MB, MDt, MC > > B
Definition: schur_matrix.h:195
void reinit(const unsigned int num_blocks, const size_type block_size=0, const bool fast=false)
void debug_level(unsigned int l)
Definition: schur_matrix.h:242
VectorMemory< BlockVector< double > > & mem
Definition: schur_matrix.h:207
SchurMatrix & operator=(const SchurMatrix< MA_inverse, MB, MDt, MC > &)
const SmartPointer< const MA_inverse, SchurMatrix< MA_inverse, MB, MDt, MC > > Ainv
Definition: schur_matrix.h:191
void sadd(const value_type s, const BlockVectorBase &V)
#define Assert(cond, exc)
Definition: exceptions.h:299
BlockType & block(const unsigned int i)
const SmartPointer< const MDt, SchurMatrix< MA_inverse, MB, MDt, MC > > Dt
Definition: schur_matrix.h:199
void push(const std::string &text)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void scale(const value_type factor) DEAL_II_DEPRECATED
const SmartPointer< const MC, SchurMatrix< MA_inverse, MB, MDt, MC > > C
Definition: schur_matrix.h:203