Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_matrix_array.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/lac/block_matrix_array.h>
18 #include <deal.II/lac/block_vector.h>
19 #include <deal.II/lac/petsc_block_vector.h>
20 #include <deal.II/lac/trilinos_parallel_block_vector.h>
21 #include <deal.II/lac/vector.h>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 
26 template <typename number, typename BlockVectorType>
28  : row(e.row)
29  , col(e.col)
30  , prefix(e.prefix)
31  , transpose(e.transpose)
32  , matrix(e.matrix)
33 {
34  Entry &e2 = const_cast<Entry &>(e);
35  e2.matrix = nullptr;
36 }
37 
38 
39 
40 template <typename number, typename BlockVectorType>
42 {
43  if (matrix)
44  delete matrix;
45 }
46 
47 
48 template <typename number, typename BlockVectorType>
50  : block_rows(0)
51  , block_cols(0)
52 {}
53 
54 
55 
56 template <typename number, typename BlockVectorType>
58  const unsigned int n_block_rows,
59  const unsigned int n_block_cols)
60  : block_rows(n_block_rows)
61  , block_cols(n_block_cols)
62 {}
63 
64 
65 template <typename number, typename BlockVectorType>
66 void
68  const unsigned int n_block_rows,
69  const unsigned int n_block_cols)
70 {
73 }
74 
75 
76 
77 template <typename number, typename BlockVectorType>
78 void
80  const unsigned int n_block_rows,
81  const unsigned int n_block_cols)
82 {
83  clear();
86 }
87 
88 
89 
90 template <typename number, typename BlockVectorType>
91 void
93 {
94  entries.clear();
95 }
96 
97 
98 template <typename number, typename BlockVectorType>
99 void
101  BlockVectorType & dst,
102  const BlockVectorType &src) const
103 {
105  Assert(dst.n_blocks() == block_rows,
106  ExcDimensionMismatch(dst.n_blocks(), block_rows));
107  Assert(src.n_blocks() == block_cols,
108  ExcDimensionMismatch(src.n_blocks(), block_cols));
109 
111  mem);
112  typename BlockVectorType::BlockType &aux = *p_aux;
113 
114  typename std::vector<Entry>::const_iterator m = entries.begin();
115  typename std::vector<Entry>::const_iterator end = entries.end();
116 
117  for (; m != end; ++m)
118  {
119  aux.reinit(dst.block(m->row));
120  if (m->transpose)
121  m->matrix->Tvmult(aux, src.block(m->col));
122  else
123  m->matrix->vmult(aux, src.block(m->col));
124  dst.block(m->row).add(m->prefix, aux);
125  }
126 }
127 
128 
129 
130 template <typename number, typename BlockVectorType>
131 void
133  BlockVectorType & dst,
134  const BlockVectorType &src) const
135 {
136  dst = 0.;
137  vmult_add(dst, src);
138 }
139 
140 
141 
142 template <typename number, typename BlockVectorType>
143 void
145  BlockVectorType & dst,
146  const BlockVectorType &src) const
147 {
149  Assert(dst.n_blocks() == block_cols,
150  ExcDimensionMismatch(dst.n_blocks(), block_cols));
151  Assert(src.n_blocks() == block_rows,
152  ExcDimensionMismatch(src.n_blocks(), block_rows));
153 
154  typename std::vector<Entry>::const_iterator m = entries.begin();
155  typename std::vector<Entry>::const_iterator end = entries.end();
156 
158  mem);
159  typename BlockVectorType::BlockType &aux = *p_aux;
160 
161  for (; m != end; ++m)
162  {
163  aux.reinit(dst.block(m->col));
164  if (m->transpose)
165  m->matrix->vmult(aux, src.block(m->row));
166  else
167  m->matrix->Tvmult(aux, src.block(m->row));
168  dst.block(m->col).add(m->prefix, aux);
169  }
170 }
171 
172 
173 
174 template <typename number, typename BlockVectorType>
175 void
177  BlockVectorType & dst,
178  const BlockVectorType &src) const
179 {
180  dst = 0.;
181  Tvmult_add(dst, src);
182 }
183 
184 
185 
186 template <typename number, typename BlockVectorType>
187 number
189  const BlockVectorType &u,
190  const BlockVectorType &v) const
191 {
193  Assert(u.n_blocks() == block_rows,
194  ExcDimensionMismatch(u.n_blocks(), block_rows));
195  Assert(v.n_blocks() == block_cols,
196  ExcDimensionMismatch(v.n_blocks(), block_cols));
197 
199  mem);
200  typename BlockVectorType::BlockType &aux = *p_aux;
201 
202  typename std::vector<Entry>::const_iterator m;
203  typename std::vector<Entry>::const_iterator end = entries.end();
204 
205  number result = 0.;
206 
207  for (unsigned int i = 0; i < block_rows; ++i)
208  {
209  aux.reinit(u.block(i));
210  for (m = entries.begin(); m != end; ++m)
211  {
212  if (m->row != i)
213  continue;
214  if (m->transpose)
215  m->matrix->Tvmult_add(aux, v.block(m->col));
216  else
217  m->matrix->vmult(aux, v.block(m->col));
218  }
219  result += u.block(i) * aux;
220  }
221 
222  return result;
223 }
224 
225 
226 
227 template <typename number, typename BlockVectorType>
228 number
230  const BlockVectorType &u) const
231 {
232  return matrix_scalar_product(u, u);
233 }
234 
235 
236 
237 template <typename number, typename BlockVectorType>
238 unsigned int
240 {
241  return block_rows;
242 }
243 
244 
245 
246 template <typename number, typename BlockVectorType>
247 unsigned int
249 {
250  return block_cols;
251 }
252 
253 
254 
255 //---------------------------------------------------------------------------
256 
257 template <typename number, typename BlockVectorType>
259  : BlockMatrixArray<number, BlockVectorType>()
260  , backward(false)
261 {}
262 
263 
264 template <typename number, typename BlockVectorType>
266  const unsigned int block_rows)
267  : BlockMatrixArray<number, BlockVectorType>(block_rows, block_rows)
268  , backward(false)
269 {}
270 
271 
272 template <typename number, typename BlockVectorType>
273 void
275 {
277 }
278 
279 
280 template <typename number, typename BlockVectorType>
281 void
283  BlockVectorType &dst,
284  size_type row_num) const
285 {
287  typename std::vector<
288  typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator
289  m = this->entries.begin();
290  typename std::vector<
291  typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator
292  end = this->entries.end();
293  std::vector<typename std::vector<
294  typename BlockMatrixArray<number, BlockVectorType>::Entry>::const_iterator>
295  diagonals;
296 
298  mem);
299  typename BlockVectorType::BlockType &aux = *p_aux;
300 
301  aux.reinit(dst.block(row_num), true);
302 
303  // Loop over all entries, since
304  // they are not ordered by rows.
305  for (; m != end; ++m)
306  {
307  const size_type i = m->row;
308  // Ignore everything not in
309  // this row
310  if (i != row_num)
311  continue;
312  const size_type j = m->col;
313  // Only use the lower (upper)
314  // triangle for forward
315  // (backward) substitution
316  if (((j > i) && !backward) || ((j < i) && backward))
317  continue;
318  if (j == i)
319  {
320  diagonals.push_back(m);
321  }
322  else
323  {
324  if (m->transpose)
325  m->matrix->Tvmult(aux, dst.block(j));
326  else
327  m->matrix->vmult(aux, dst.block(j));
328  dst.block(i).add(-1.0 * m->prefix, aux);
329  }
330  }
331  Assert(diagonals.size() != 0, ExcNoDiagonal(row_num));
332 
333  // Inverting the diagonal block is
334  // simple, if there is only one
335  // matrix
336  if (diagonals.size() == 1)
337  {
338  if (diagonals[0]->transpose)
339  diagonals[0]->matrix->Tvmult(aux, dst.block(row_num));
340  else
341  diagonals[0]->matrix->vmult(aux, dst.block(row_num));
342  dst.block(row_num).equ(diagonals[0]->prefix, aux);
343  }
344  else
345  {
346  aux = 0.;
347  for (size_type i = 0; i < diagonals.size(); ++i)
348  {
349  m = diagonals[i];
350  // First, divide by the current
351  // factor, such that we can
352  // multiply by it later.
353  aux /= m->prefix;
354  if (m->transpose)
355  m->matrix->Tvmult_add(aux, dst.block(row_num));
356  else
357  m->matrix->vmult_add(aux, dst.block(row_num));
358  aux *= m->prefix;
359  }
360  dst.block(row_num) = aux;
361  }
362 }
363 
364 
365 
366 template <typename number, typename BlockVectorType>
367 void
369  BlockVectorType & dst,
370  const BlockVectorType &src) const
371 {
372  Assert(dst.n_blocks() == n_block_rows(),
373  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
374  Assert(src.n_blocks() == n_block_cols(),
375  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
376 
377  BlockVectorType aux;
378  aux.reinit(dst);
379  vmult(aux, src);
380  dst += aux;
381 }
382 
383 
384 
385 template <typename number, typename BlockVectorType>
386 void
388  BlockVectorType & dst,
389  const BlockVectorType &src) const
390 {
391  Assert(dst.n_blocks() == n_block_rows(),
392  ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
393  Assert(src.n_blocks() == n_block_cols(),
394  ExcDimensionMismatch(src.n_blocks(), n_block_cols()));
395 
396  dst.equ(1., src);
397 
398  if (backward)
399  {
400  for (unsigned int i = n_block_rows(); i > 0;)
401  do_row(dst, --i);
402  }
403  else
404  {
405  for (unsigned int i = 0; i < n_block_rows(); ++i)
406  do_row(dst, i);
407  }
408 }
409 
410 template <typename number, typename BlockVectorType>
411 void
413  BlockVectorType &,
414  const BlockVectorType &) const
415 {
416  Assert(false, ExcNotImplemented());
417 }
418 
419 
420 template <typename number, typename BlockVectorType>
421 void
423  BlockVectorType &,
424  const BlockVectorType &) const
425 {
426  Assert(false, ExcNotImplemented());
427 }
428 
429 template class BlockMatrixArray<float>;
430 template class BlockMatrixArray<double>;
431 template class BlockTrianglePrecondition<float>;
432 template class BlockTrianglePrecondition<double>;
433 
434 #ifdef DEAL_II_WITH_TRILINOS
437 template class BlockTrianglePrecondition<float,
439 template class BlockTrianglePrecondition<double,
440  TrilinosWrappers::MPI::BlockVector>;
441 #endif
442 
443 #ifdef DEAL_II_WITH_PETSC
444 # ifdef PETSC_USE_COMPLEX
448  PETScWrappers::MPI::BlockVector>;
450  PETScWrappers::MPI::BlockVector>;
452  PETScWrappers::MPI::BlockVector>;
453 # else
456 template class BlockTrianglePrecondition<float,
457  PETScWrappers::MPI::BlockVector>;
458 template class BlockTrianglePrecondition<double,
459  PETScWrappers::MPI::BlockVector>;
460 # endif
461 #endif
462 
463 DEAL_II_NAMESPACE_CLOSE
unsigned int block_cols
std::vector< Entry > entries
number matrix_norm_square(const BlockVectorType &u) const
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
PointerMatrixBase< typename BlockVectorType::BlockType > * matrix
void Tvmult(BlockVectorType &dst, const BlockVectorType &src) const
void vmult(BlockVectorType &dst, const BlockVectorType &src) const
types::global_dof_index size_type
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
unsigned int n_block_cols() const
void initialize(const unsigned int n_block_rows, const unsigned int n_block_cols)
#define Assert(cond, exc)
Definition: exceptions.h:1407
void vmult_add(BlockVectorType &dst, const BlockVectorType &src) const
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
unsigned int n_block_rows() const
void Tvmult(BlockVectorType &dst, const BlockVectorType &src) const
unsigned int block_rows
Entry(const MatrixType &matrix, size_type row, size_type col, number prefix, bool transpose)
void reinit(const unsigned int n_block_rows, const unsigned int n_block_cols)
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
static ::ExceptionBase & ExcNoDiagonal(size_type arg1)
void Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const
static ::ExceptionBase & ExcNotImplemented()
void reinit(const unsigned int n_block_rows)
void vmult(BlockVectorType &dst, const BlockVectorType &src) const
void do_row(BlockVectorType &dst, size_type row_num) const
number matrix_scalar_product(const BlockVectorType &u, const BlockVectorType &v) const
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)