Reference documentation for deal.II version Git 437f94a 2018-06-23 09:16:19 -0400
block_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2017 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/base/table.h>
23 
24 #include <deal.II/lac/block_matrix_base.h>
25 #include <deal.II/lac/block_sparsity_pattern.h>
26 #include <deal.II/lac/block_vector.h>
27 #include <deal.II/lac/exceptions.h>
28 #include <deal.II/lac/sparse_matrix.h>
29 
30 #include <cmath>
31 
32 DEAL_II_NAMESPACE_OPEN
33 
34 
51 template <typename number>
52 class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix<number>>
53 {
54 public:
59 
63  typedef typename BaseClass::BlockType BlockType;
64 
69  typedef typename BaseClass::pointer pointer;
70  typedef typename BaseClass::const_pointer const_pointer;
71  typedef typename BaseClass::reference reference;
72  typedef typename BaseClass::const_reference const_reference;
73  typedef typename BaseClass::size_type size_type;
74  typedef typename BaseClass::iterator iterator;
76 
92  BlockSparseMatrix() = default;
93 
106  BlockSparseMatrix(const BlockSparsityPattern &sparsity);
107 
111  virtual ~BlockSparseMatrix() override;
112 
113 
114 
120  operator=(const BlockSparseMatrix &);
121 
131  operator=(const double d);
132 
141  void
142  clear();
143 
158  virtual void
159  reinit(const BlockSparsityPattern &sparsity);
161 
170  bool
171  empty() const;
172 
176  size_type
177  get_row_length(const size_type row) const;
178 
184  size_type
185  n_nonzero_elements() const;
186 
192  size_type
193  n_actually_nonzero_elements(const double threshold = 0.0) const;
194 
203  const BlockSparsityPattern &
204  get_sparsity_pattern() const;
205 
210  std::size_t
211  memory_consumption() const;
213 
222  template <typename block_number>
223  void
225  const BlockVector<block_number> &src) const;
226 
231  template <typename block_number, typename nonblock_number>
232  void
234  const Vector<nonblock_number> &src) const;
235 
240  template <typename block_number, typename nonblock_number>
241  void
242  vmult(Vector<nonblock_number> & dst,
243  const BlockVector<block_number> &src) const;
244 
249  template <typename nonblock_number>
250  void
251  vmult(Vector<nonblock_number> &dst, const Vector<nonblock_number> &src) const;
252 
258  template <typename block_number>
259  void
261  const BlockVector<block_number> &src) const;
262 
267  template <typename block_number, typename nonblock_number>
268  void
270  const Vector<nonblock_number> &src) const;
271 
276  template <typename block_number, typename nonblock_number>
277  void
278  Tvmult(Vector<nonblock_number> & dst,
279  const BlockVector<block_number> &src) const;
280 
285  template <typename nonblock_number>
286  void
287  Tvmult(Vector<nonblock_number> & dst,
288  const Vector<nonblock_number> &src) const;
290 
302  template <class BlockVectorType>
303  void
304  precondition_Jacobi(BlockVectorType & dst,
305  const BlockVectorType &src,
306  const number omega = 1.) const;
307 
313  template <typename number2>
314  void
315  precondition_Jacobi(Vector<number2> & dst,
316  const Vector<number2> &src,
317  const number omega = 1.) const;
319 
344  void
345  print_formatted(std::ostream & out,
346  const unsigned int precision = 3,
347  const bool scientific = true,
348  const unsigned int width = 0,
349  const char * zero_string = " ",
350  const double denominator = 1.) const;
352 
362 
363 private:
371 };
372 
373 
374 
376 /* ------------------------- Template functions ---------------------- */
377 
378 
379 
380 template <typename number>
383 {
385 
386  for (size_type r = 0; r < this->n_block_rows(); ++r)
387  for (size_type c = 0; c < this->n_block_cols(); ++c)
388  this->block(r, c) = d;
389 
390  return *this;
391 }
392 
393 
394 
395 template <typename number>
396 template <typename block_number>
397 inline void
399  const BlockVector<block_number> &src) const
400 {
401  BaseClass::vmult_block_block(dst, src);
402 }
403 
404 
405 
406 template <typename number>
407 template <typename block_number, typename nonblock_number>
408 inline void
410  const Vector<nonblock_number> &src) const
411 {
412  BaseClass::vmult_block_nonblock(dst, src);
413 }
414 
415 
416 
417 template <typename number>
418 template <typename block_number, typename nonblock_number>
419 inline void
420 BlockSparseMatrix<number>::vmult(Vector<nonblock_number> & dst,
421  const BlockVector<block_number> &src) const
422 {
423  BaseClass::vmult_nonblock_block(dst, src);
424 }
425 
426 
427 
428 template <typename number>
429 template <typename nonblock_number>
430 inline void
431 BlockSparseMatrix<number>::vmult(Vector<nonblock_number> & dst,
432  const Vector<nonblock_number> &src) const
433 {
434  BaseClass::vmult_nonblock_nonblock(dst, src);
435 }
436 
437 
438 
439 template <typename number>
440 template <typename block_number>
441 inline void
443  const BlockVector<block_number> &src) const
444 {
445  BaseClass::Tvmult_block_block(dst, src);
446 }
447 
448 
449 
450 template <typename number>
451 template <typename block_number, typename nonblock_number>
452 inline void
454  const Vector<nonblock_number> &src) const
455 {
456  BaseClass::Tvmult_block_nonblock(dst, src);
457 }
458 
459 
460 
461 template <typename number>
462 template <typename block_number, typename nonblock_number>
463 inline void
464 BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> & dst,
465  const BlockVector<block_number> &src) const
466 {
467  BaseClass::Tvmult_nonblock_block(dst, src);
468 }
469 
470 
471 
472 template <typename number>
473 template <typename nonblock_number>
474 inline void
475 BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> & dst,
476  const Vector<nonblock_number> &src) const
477 {
478  BaseClass::Tvmult_nonblock_nonblock(dst, src);
479 }
480 
481 
482 
483 template <typename number>
484 template <class BlockVectorType>
485 inline void
487  const BlockVectorType &src,
488  const number omega) const
489 {
490  Assert(this->n_block_rows() == this->n_block_cols(), ExcNotQuadratic());
491  Assert(dst.n_blocks() == this->n_block_rows(),
492  ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows()));
493  Assert(src.n_blocks() == this->n_block_cols(),
494  ExcDimensionMismatch(src.n_blocks(), this->n_block_cols()));
495 
496  // do a diagonal preconditioning. uses only
497  // the diagonal blocks of the matrix
498  for (size_type i = 0; i < this->n_block_rows(); ++i)
499  this->block(i, i).precondition_Jacobi(dst.block(i), src.block(i), omega);
500 }
501 
502 
503 
504 template <typename number>
505 template <typename number2>
506 inline void
508  const Vector<number2> &src,
509  const number omega) const
510 {
511  // check number of blocks. the sizes of the
512  // single block is checked in the function
513  // we call
514  Assert(this->n_block_cols() == 1,
515  ExcMessage("This function only works if the matrix has "
516  "a single block"));
517  Assert(this->n_block_rows() == 1,
518  ExcMessage("This function only works if the matrix has "
519  "a single block"));
520 
521  // do a diagonal preconditioning. uses only
522  // the diagonal blocks of the matrix
523  this->block(0, 0).precondition_Jacobi(dst, src, omega);
524 }
525 
526 
527 DEAL_II_NAMESPACE_CLOSE
528 
529 #endif // dealii_block_sparse_matrix_h
SmartPointer< const BlockSparsityPattern, BlockSparseMatrix< number > > sparsity_pattern
virtual ~BlockSparseMatrix() override
static::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
std::size_t memory_consumption() const
void precondition_Jacobi(BlockVectorType &dst, const BlockVectorType &src, const number omega=1.) const
size_type n_actually_nonzero_elements(const double threshold=0.0) const
BaseClass::value_type value_type
BlockSparseMatrix & operator=(const BlockSparseMatrix &)
size_type get_row_length(const size_type row) const
void Tvmult(BlockVector< block_number > &dst, const BlockVector< block_number > &src) const
static::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1227
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:385
void vmult(BlockVector< block_number > &dst, const BlockVector< block_number > &src) 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
size_type n_nonzero_elements() const
BlockMatrixBase< SparseMatrix< number > > BaseClass
static::ExceptionBase & ExcNotQuadratic()
const BlockSparsityPattern & get_sparsity_pattern() const
BlockSparseMatrix()=default
BaseClass::BlockType BlockType
bool empty() const
static::ExceptionBase & ExcBlockDimensionMismatch()
virtual void reinit(const BlockSparsityPattern &sparsity)