Reference documentation for deal.II version 8.1.0
relaxation_block.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: relaxation_block.templates.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 1999 - 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__relaxation_block_templates_h
18 #define __deal2__relaxation_block_templates_h
19 
20 #include <deal.II/lac/relaxation_block.h>
21 #include <deal.II/lac/full_matrix.h>
22 #include <deal.II/lac/vector_memory.h>
23 
25 
26 template <class MATRIX, typename inverse_type>
27 inline
29  const double relaxation,
30  const bool invert_diagonal,
31  const bool same_diagonal)
32  :
33  relaxation(relaxation),
34  invert_diagonal(invert_diagonal),
35  same_diagonal(same_diagonal),
36  inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
37  threshold(0.)
38 {}
39 
40 
41 template <class MATRIX, typename inverse_type>
42 inline
43 std::size_t
45 {
46  std::size_t result = sizeof(*this)
47  + - sizeof(block_list) + block_list.memory_consumption();
48  for (unsigned int i=0; i<order.size(); ++i)
49  result += MemoryConsumption::memory_consumption(order[i]);
50  return result;
51 }
52 
53 
54 template <class MATRIX, typename inverse_type>
55 inline
57  const BlockList &bl,
58  const double relaxation,
59  const bool invert_diagonal,
60  const bool same_diagonal)
61  :
62  relaxation(relaxation),
63  invert_diagonal(invert_diagonal),
64  same_diagonal(same_diagonal),
66  threshold(0.)
67 {
69 }
70 
71 
72 template <class MATRIX, typename inverse_type>
73 inline
74 void
76  const MATRIX &M,
77  const AdditionalData &parameters)
78 {
80 
81  clear();
82 // Assert (M.m() == M.n(), ExcNotQuadratic());
83  A = &M;
84  additional_data = &parameters;
85  this->inversion = parameters.inversion;
86 
87  this->reinit(additional_data->block_list.n_rows(), 0, additional_data->same_diagonal,
88  additional_data->inversion);
89 
90  if (additional_data->invert_diagonal)
92 }
93 
94 
95 template <class MATRIX, typename inverse_type>
96 inline
97 void
99 {
100  A = 0;
102 }
103 
104 
105 template <class MATRIX, typename inverse_type>
106 inline
107 void
109 {
110  const MATRIX &M=*A;
112 
113  if (this->same_diagonal())
114  {
115  Assert(false, ExcNotImplemented());
116  }
117  else
118  {
119  for (size_type block=0; block<additional_data->block_list.n_rows(); ++block)
120  {
121  const size_type bs = additional_data->block_list.row_length(block);
122  M_cell.reinit(bs, bs);
123 
124  // Copy rows for this block
125  // into the matrix for the
126  // diagonal block
128  = additional_data->block_list.begin(block);
129  for (size_type row_cell=0; row_cell<bs; ++row_cell, ++row)
130  {
131 //TODO:[GK] Optimize here
132  for (typename MATRIX::const_iterator entry = M.begin(row->column());
133  entry != M.end(row->column()); ++entry)
134  {
135  const size_type column = entry->column();
136  const size_type col_cell = additional_data->block_list.row_position(block, column);
137  if (col_cell != numbers::invalid_size_type)
138  M_cell(row_cell, col_cell) = entry->value();
139  }
140  }
141  // Now M_cell contains the
142  // diagonal block. Now
143  // store it and its
144  // inverse, if so requested.
145  if (this->store_diagonals())
146  {
147  this->diagonal(block).reinit(bs, bs);
148  this->diagonal(block) = M_cell;
149  }
150  switch (this->inversion)
151  {
153  this->inverse(block).reinit(bs, bs);
154  this->inverse(block).invert(M_cell);
155  break;
157  this->inverse_householder(block).initialize(M_cell);
158  break;
160  this->inverse_svd(block).reinit(bs, bs);
161  this->inverse_svd(block) = M_cell;
162  this->inverse_svd(block).compute_inverse_svd(additional_data->threshold);
163  break;
164  default:
165  Assert(false, ExcNotImplemented());
166  }
167  }
168  }
169  this->inverses_computed(true);
170 }
171 
172 
173 template <class MATRIX, typename inverse_type>
174 template <typename number2>
175 inline
176 void
178  Vector<number2> &dst,
179  const Vector<number2> &prev,
180  const Vector<number2> &src,
181  const bool backward) const
182 {
183  Assert (additional_data->invert_diagonal, ExcNotImplemented());
184 
185  const MATRIX &M=*this->A;
186  Vector<number2> b_cell, x_cell;
187 
188  const bool permutation_empty = additional_data->order.size() == 0;
189  const unsigned int n_permutations = (permutation_empty)
190  ? 1U : additional_data->order.size();
191  const size_type n_blocks = additional_data->block_list.n_rows();
192 
193  if (!permutation_empty)
194  for (unsigned int i=0; i<additional_data->order.size(); ++i)
195  AssertDimension(additional_data->order[i].size(), this->size());
196 
197  for (unsigned int perm=0; perm<n_permutations; ++perm)
198  {
199  for (unsigned int bi=0; bi<n_blocks; ++bi)
200  {
201  const unsigned int raw_block = backward ? (n_blocks - bi - 1) : bi;
202  const unsigned int block = permutation_empty
203  ? raw_block
204  : (backward
205  ? (additional_data->order[n_permutations-1-perm][raw_block])
206  : (additional_data->order[perm][raw_block]));
207 
208  const size_type bs = additional_data->block_list.row_length(block);
209 
210  b_cell.reinit(bs);
211  x_cell.reinit(bs);
212  // Collect off-diagonal parts
213  SparsityPattern::iterator row = additional_data->block_list.begin(block);
214  for (size_type row_cell=0; row_cell<bs; ++row_cell, ++row)
215  {
216  b_cell(row_cell) = src(row->column());
217  for (typename MATRIX::const_iterator entry = M.begin(row->column());
218  entry != M.end(row->column()); ++entry)
219  b_cell(row_cell) -= entry->value() * prev(entry->column());
220  }
221  // Apply inverse diagonal
222  this->inverse_vmult(block, x_cell, b_cell);
223 #ifdef DEBUG
224  for (unsigned int i=0; i<x_cell.size(); ++i)
225  {
227  }
228 #endif
229  // Store in result vector
230  row=additional_data->block_list.begin(block);
231  for (size_type row_cell=0; row_cell<bs; ++row_cell, ++row)
232  dst(row->column()) = prev(row->column()) + additional_data->relaxation * x_cell(row_cell);
233  }
234  }
235 }
236 
237 
238 //----------------------------------------------------------------------//
239 
240 template <class MATRIX, typename inverse_type>
241 template <typename number2>
243  Vector<number2> &dst,
244  const Vector<number2> &src) const
245 {
247  typename VectorMemory<Vector<number2> >::Pointer aux = mem;
248  aux->reinit(dst, false);
249  *aux = dst;
250  this->do_step(dst, *aux, src, false);
251 }
252 
253 
254 template <class MATRIX, typename inverse_type>
255 template <typename number2>
257  Vector<number2> &dst,
258  const Vector<number2> &src) const
259 {
261  typename VectorMemory<Vector<number2> >::Pointer aux = mem;
262  aux->reinit(dst, false);
263  *aux = dst;
264  this->do_step(dst, *aux, src, true);
265 }
266 
267 
268 //----------------------------------------------------------------------//
269 
270 template <class MATRIX, typename inverse_type>
271 template <typename number2>
273  Vector<number2> &dst,
274  const Vector<number2> &src) const
275 {
276  this->do_step(dst, dst, src, false);
277 }
278 
279 
280 template <class MATRIX, typename inverse_type>
281 template <typename number2>
283  Vector<number2> &dst,
284  const Vector<number2> &src) const
285 {
286  this->do_step(dst, dst, src, true);
287 }
288 
289 
290 //----------------------------------------------------------------------//
291 
292 template <class MATRIX, typename inverse_type>
293 template <typename number2>
295  Vector<number2> &dst,
296  const Vector<number2> &src) const
297 {
298  this->do_step(dst, dst, src, false);
299  this->do_step(dst, dst, src, true);
300 }
301 
302 
303 template <class MATRIX, typename inverse_type>
304 template <typename number2>
306  Vector<number2> &dst,
307  const Vector<number2> &src) const
308 {
309  this->do_step(dst, dst, src, true);
310  this->do_step(dst, dst, src, false);
311 }
312 
313 
314 
315 DEAL_II_NAMESPACE_CLOSE
316 
317 
318 #endif
const types::global_dof_index invalid_size_type
Definition: types.h:211
void reinit(const TableIndices< N > &new_size, const bool fast=false)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
SmartPointer< const MATRIX, RelaxationBlock< MATRIX, inverse_type > > A
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void invert(const FullMatrix< number2 > &M)
bool is_finite(const double x)
FullMatrix< inverse_type > & inverse(size_type i)
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
SmartPointer< const AdditionalData, RelaxationBlock< MATRIX, inverse_type > > additional_data
std::size_t size() const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void reinit(unsigned int nblocks, size_type blocksize, bool compress, Inversion method=gauss_jordan)
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void create_sparsity_pattern(SparsityPattern &sparsity, size_type n) const
Definition: block_list.h:302
void compute_inverse_svd(const double threshold=0.)
::ExceptionBase & ExcNumberNotFinite()
void initialize(const MATRIX &A, const AdditionalData &parameters)
void do_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool backward) const
PreconditionBlockBase< inverse_type >::Inversion inversion
void reinit(const unsigned int size1, const unsigned int size2, const bool fast=false)
AdditionalData(const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
LAPACKFullMatrix< inverse_type > & inverse_svd(size_type i)
Householder< inverse_type > & inverse_householder(size_type i)
::ExceptionBase & ExcNotImplemented()
void initialize(const FullMatrix< number2 > &)
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
FullMatrix< inverse_type > & diagonal(size_type i)
types::global_dof_index size_type
void inverse_vmult(size_type i, Vector< number2 > &dst, const Vector< number2 > &src) const