Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 #ifndef dealii_block_sparse_matrix_h
17 #define dealii_block_sparse_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/lac/block_matrix_base.h>
23 #include <deal.II/lac/block_sparsity_pattern.h>
24 #include <deal.II/lac/block_vector.h>
25 #include <deal.II/lac/exceptions.h>
26 #include <deal.II/lac/sparse_matrix.h>
27 
28 #include <cmath>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
49 template <typename number>
50 class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix<number>>
51 {
52 public:
57 
61  using BlockType = typename BaseClass::BlockType;
62 
67  using pointer = typename BaseClass::pointer;
68  using const_pointer = typename BaseClass::const_pointer;
69  using reference = typename BaseClass::reference;
70  using const_reference = typename BaseClass::const_reference;
71  using size_type = typename BaseClass::size_type;
72  using iterator = typename BaseClass::iterator;
73  using const_iterator = typename BaseClass::const_iterator;
74 
90  BlockSparseMatrix() = default;
91 
104  BlockSparseMatrix(const BlockSparsityPattern &sparsity);
105 
109  virtual ~BlockSparseMatrix() override;
110 
111 
112 
118  operator=(const BlockSparseMatrix &);
119 
129  operator=(const double d);
130 
139  void
140  clear();
141 
156  virtual void
157  reinit(const BlockSparsityPattern &sparsity);
159 
168  bool
169  empty() const;
170 
174  size_type
175  get_row_length(const size_type row) const;
176 
182  size_type
183  n_nonzero_elements() const;
184 
190  size_type
191  n_actually_nonzero_elements(const double threshold = 0.0) const;
192 
201  const BlockSparsityPattern &
202  get_sparsity_pattern() const;
203 
208  std::size_t
209  memory_consumption() const;
211 
220  template <typename block_number>
221  void
223  const BlockVector<block_number> &src) const;
224 
229  template <typename block_number, typename nonblock_number>
230  void
232  const Vector<nonblock_number> &src) const;
233 
238  template <typename block_number, typename nonblock_number>
239  void
240  vmult(Vector<nonblock_number> & dst,
241  const BlockVector<block_number> &src) const;
242 
247  template <typename nonblock_number>
248  void
249  vmult(Vector<nonblock_number> &dst, const Vector<nonblock_number> &src) const;
250 
256  template <typename block_number>
257  void
259  const BlockVector<block_number> &src) const;
260 
265  template <typename block_number, typename nonblock_number>
266  void
268  const Vector<nonblock_number> &src) const;
269 
274  template <typename block_number, typename nonblock_number>
275  void
276  Tvmult(Vector<nonblock_number> & dst,
277  const BlockVector<block_number> &src) const;
278 
283  template <typename nonblock_number>
284  void
285  Tvmult(Vector<nonblock_number> & dst,
286  const Vector<nonblock_number> &src) const;
288 
300  template <class BlockVectorType>
301  void
302  precondition_Jacobi(BlockVectorType & dst,
303  const BlockVectorType &src,
304  const number omega = 1.) const;
305 
311  template <typename number2>
312  void
313  precondition_Jacobi(Vector<number2> & dst,
314  const Vector<number2> &src,
315  const number omega = 1.) const;
317 
342  void
343  print_formatted(std::ostream & out,
344  const unsigned int precision = 3,
345  const bool scientific = true,
346  const unsigned int width = 0,
347  const char * zero_string = " ",
348  const double denominator = 1.) const;
350 
360 
361 private:
369 };
370 
371 
372 
374 /* ------------------------- Template functions ---------------------- */
375 
376 
377 
378 template <typename number>
381 {
383 
384  for (size_type r = 0; r < this->n_block_rows(); ++r)
385  for (size_type c = 0; c < this->n_block_cols(); ++c)
386  this->block(r, c) = d;
387 
388  return *this;
389 }
390 
391 
392 
393 template <typename number>
394 template <typename block_number>
395 inline void
397  const BlockVector<block_number> &src) const
398 {
399  BaseClass::vmult_block_block(dst, src);
400 }
401 
402 
403 
404 template <typename number>
405 template <typename block_number, typename nonblock_number>
406 inline void
408  const Vector<nonblock_number> &src) const
409 {
410  BaseClass::vmult_block_nonblock(dst, src);
411 }
412 
413 
414 
415 template <typename number>
416 template <typename block_number, typename nonblock_number>
417 inline void
418 BlockSparseMatrix<number>::vmult(Vector<nonblock_number> & dst,
419  const BlockVector<block_number> &src) const
420 {
421  BaseClass::vmult_nonblock_block(dst, src);
422 }
423 
424 
425 
426 template <typename number>
427 template <typename nonblock_number>
428 inline void
429 BlockSparseMatrix<number>::vmult(Vector<nonblock_number> & dst,
430  const Vector<nonblock_number> &src) const
431 {
432  BaseClass::vmult_nonblock_nonblock(dst, src);
433 }
434 
435 
436 
437 template <typename number>
438 template <typename block_number>
439 inline void
441  const BlockVector<block_number> &src) const
442 {
443  BaseClass::Tvmult_block_block(dst, src);
444 }
445 
446 
447 
448 template <typename number>
449 template <typename block_number, typename nonblock_number>
450 inline void
452  const Vector<nonblock_number> &src) const
453 {
454  BaseClass::Tvmult_block_nonblock(dst, src);
455 }
456 
457 
458 
459 template <typename number>
460 template <typename block_number, typename nonblock_number>
461 inline void
462 BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> & dst,
463  const BlockVector<block_number> &src) const
464 {
465  BaseClass::Tvmult_nonblock_block(dst, src);
466 }
467 
468 
469 
470 template <typename number>
471 template <typename nonblock_number>
472 inline void
473 BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> & dst,
474  const Vector<nonblock_number> &src) const
475 {
476  BaseClass::Tvmult_nonblock_nonblock(dst, src);
477 }
478 
479 
480 
481 template <typename number>
482 template <class BlockVectorType>
483 inline void
485  const BlockVectorType &src,
486  const number omega) const
487 {
488  Assert(this->n_block_rows() == this->n_block_cols(), ExcNotQuadratic());
489  Assert(dst.n_blocks() == this->n_block_rows(),
490  ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows()));
491  Assert(src.n_blocks() == this->n_block_cols(),
492  ExcDimensionMismatch(src.n_blocks(), this->n_block_cols()));
493 
494  // do a diagonal preconditioning. uses only
495  // the diagonal blocks of the matrix
496  for (size_type i = 0; i < this->n_block_rows(); ++i)
497  this->block(i, i).precondition_Jacobi(dst.block(i), src.block(i), omega);
498 }
499 
500 
501 
502 template <typename number>
503 template <typename number2>
504 inline void
506  const Vector<number2> &src,
507  const number omega) const
508 {
509  // check number of blocks. the sizes of the
510  // single block is checked in the function
511  // we call
512  Assert(this->n_block_cols() == 1,
513  ExcMessage("This function only works if the matrix has "
514  "a single block"));
515  Assert(this->n_block_rows() == 1,
516  ExcMessage("This function only works if the matrix has "
517  "a single block"));
518 
519  // do a diagonal preconditioning. uses only
520  // the diagonal blocks of the matrix
521  this->block(0, 0).precondition_Jacobi(dst, src, omega);
522 }
523 
524 
525 DEAL_II_NAMESPACE_CLOSE
526 
527 #endif // dealii_block_sparse_matrix_h
SmartPointer< const BlockSparsityPattern, BlockSparseMatrix< number > > sparsity_pattern
bool empty() const
virtual ~BlockSparseMatrix() override
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void Tvmult(BlockVector< block_number > &dst, const BlockVector< block_number > &src) const
BlockSparseMatrix & operator=(const BlockSparseMatrix &)
const BlockSparsityPattern & get_sparsity_pattern() const
typename BlockType::value_type value_type
static ::ExceptionBase & ExcMessage(std::string arg1)
typename BaseClass::BlockType BlockType
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::size_t memory_consumption() const
size_type n_actually_nonzero_elements(const double threshold=0.0) const
static ::ExceptionBase & ExcNotQuadratic()
BlockSparseMatrix()=default
typename BaseClass::value_type value_type
size_type n_nonzero_elements() const
static ::ExceptionBase & ExcBlockDimensionMismatch()
virtual void reinit(const BlockSparsityPattern &sparsity)
void vmult(BlockVector< block_number > &dst, const BlockVector< block_number > &src) const
size_type get_row_length(const size_type row) const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
void precondition_Jacobi(BlockVectorType &dst, const BlockVectorType &src, const number omega=1.) const