Reference documentation for deal.II version 9.0.0
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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/table.h>
22 #include <deal.II/lac/block_matrix_base.h>
23 #include <deal.II/lac/block_vector.h>
24 #include <deal.II/lac/sparse_matrix.h>
25 #include <deal.II/lac/block_sparsity_pattern.h>
26 #include <deal.II/lac/exceptions.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  typedef typename BaseClass::BlockType BlockType;
62 
67  typedef typename BaseClass::pointer pointer;
68  typedef typename BaseClass::const_pointer const_pointer;
69  typedef typename BaseClass::reference reference;
70  typedef typename BaseClass::const_reference const_reference;
71  typedef typename BaseClass::size_type size_type;
72  typedef typename BaseClass::iterator iterator;
74 
90  BlockSparseMatrix () = default;
91 
104  BlockSparseMatrix (const BlockSparsityPattern &sparsity);
105 
109  virtual ~BlockSparseMatrix ();
110 
111 
112 
119 
129  operator = (const double d);
130 
139  void clear ();
140 
155  virtual void reinit (const BlockSparsityPattern &sparsity);
157 
166  bool empty () const;
167 
171  size_type get_row_length (const size_type row) const;
172 
178  size_type n_nonzero_elements () const;
179 
185  size_type n_actually_nonzero_elements (const double threshold = 0.0) const;
186 
195  const BlockSparsityPattern &
196  get_sparsity_pattern () const;
197 
202  std::size_t memory_consumption () const;
204 
213  template <typename block_number>
214  void vmult (BlockVector<block_number> &dst,
215  const BlockVector<block_number> &src) const;
216 
221  template <typename block_number,
222  typename nonblock_number>
223  void vmult (BlockVector<block_number> &dst,
224  const Vector<nonblock_number> &src) const;
225 
230  template <typename block_number,
231  typename nonblock_number>
232  void vmult (Vector<nonblock_number> &dst,
233  const BlockVector<block_number> &src) const;
234 
239  template <typename nonblock_number>
240  void vmult (Vector<nonblock_number> &dst,
241  const Vector<nonblock_number> &src) const;
242 
248  template <typename block_number>
250  const BlockVector<block_number> &src) const;
251 
256  template <typename block_number,
257  typename nonblock_number>
259  const Vector<nonblock_number> &src) const;
260 
265  template <typename block_number,
266  typename nonblock_number>
267  void Tvmult (Vector<nonblock_number> &dst,
268  const BlockVector<block_number> &src) const;
269 
274  template <typename nonblock_number>
275  void Tvmult (Vector<nonblock_number> &dst,
276  const Vector<nonblock_number> &src) const;
278 
290  template <class BlockVectorType>
291  void precondition_Jacobi (BlockVectorType &dst,
292  const BlockVectorType &src,
293  const number omega = 1.) const;
294 
300  template <typename number2>
301  void precondition_Jacobi (Vector<number2> &dst,
302  const Vector<number2> &src,
303  const number omega = 1.) const;
305 
330  void print_formatted (std::ostream &out,
331  const unsigned int precision = 3,
332  const bool scientific = true,
333  const unsigned int width = 0,
334  const char *zero_string = " ",
335  const double denominator = 1.) const;
337 
347 
348 private:
355 };
356 
357 
358 
360 /* ------------------------- Template functions ---------------------- */
361 
362 
363 
364 template <typename number>
365 inline
368 {
370 
371  for (size_type r=0; r<this->n_block_rows(); ++r)
372  for (size_type c=0; c<this->n_block_cols(); ++c)
373  this->block(r,c) = d;
374 
375  return *this;
376 }
377 
378 
379 
380 template <typename number>
381 template <typename block_number>
382 inline
383 void
385  const BlockVector<block_number> &src) const
386 {
387  BaseClass::vmult_block_block (dst, src);
388 }
389 
390 
391 
392 template <typename number>
393 template <typename block_number,
394  typename nonblock_number>
395 inline
396 void
398  const Vector<nonblock_number> &src) const
399 {
400  BaseClass::vmult_block_nonblock (dst, src);
401 }
402 
403 
404 
405 template <typename number>
406 template <typename block_number,
407  typename nonblock_number>
408 inline
409 void
410 BlockSparseMatrix<number>::vmult (Vector<nonblock_number> &dst,
411  const BlockVector<block_number> &src) const
412 {
413  BaseClass::vmult_nonblock_block (dst, src);
414 }
415 
416 
417 
418 template <typename number>
419 template <typename nonblock_number>
420 inline
421 void
422 BlockSparseMatrix<number>::vmult (Vector<nonblock_number> &dst,
423  const Vector<nonblock_number> &src) const
424 {
425  BaseClass::vmult_nonblock_nonblock (dst, src);
426 }
427 
428 
429 
430 template <typename number>
431 template <typename block_number>
432 inline
433 void
435  const BlockVector<block_number> &src) const
436 {
437  BaseClass::Tvmult_block_block (dst, src);
438 }
439 
440 
441 
442 template <typename number>
443 template <typename block_number,
444  typename nonblock_number>
445 inline
446 void
448  const Vector<nonblock_number> &src) const
449 {
450  BaseClass::Tvmult_block_nonblock (dst, src);
451 }
452 
453 
454 
455 template <typename number>
456 template <typename block_number,
457  typename nonblock_number>
458 inline
459 void
460 BlockSparseMatrix<number>::Tvmult (Vector<nonblock_number> &dst,
461  const BlockVector<block_number> &src) const
462 {
463  BaseClass::Tvmult_nonblock_block (dst, src);
464 }
465 
466 
467 
468 template <typename number>
469 template <typename nonblock_number>
470 inline
471 void
472 BlockSparseMatrix<number>::Tvmult (Vector<nonblock_number> &dst,
473  const Vector<nonblock_number> &src) const
474 {
475  BaseClass::Tvmult_nonblock_nonblock (dst, src);
476 }
477 
478 
479 
480 template <typename number>
481 template <class BlockVectorType>
482 inline
483 void
485 precondition_Jacobi (BlockVectorType &dst,
486  const BlockVectorType &src,
487  const number omega) const
488 {
489  Assert (this->n_block_rows() == this->n_block_cols(), ExcNotQuadratic());
490  Assert (dst.n_blocks() == this->n_block_rows(),
491  ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows()));
492  Assert (src.n_blocks() == this->n_block_cols(),
493  ExcDimensionMismatch(src.n_blocks(), this->n_block_cols()));
494 
495  // do a diagonal preconditioning. uses only
496  // the diagonal blocks of the matrix
497  for (size_type i=0; i<this->n_block_rows(); ++i)
498  this->block(i,i).precondition_Jacobi (dst.block(i),
499  src.block(i),
500  omega);
501 }
502 
503 
504 
505 template <typename number>
506 template <typename number2>
507 inline
508 void
510 precondition_Jacobi (Vector<number2> &dst,
511  const Vector<number2> &src,
512  const number omega) const
513 {
514  // check number of blocks. the sizes of the
515  // single block is checked in the function
516  // we call
517  Assert (this->n_block_cols() == 1,
518  ExcMessage ("This function only works if the matrix has "
519  "a single block"));
520  Assert (this->n_block_rows() == 1,
521  ExcMessage ("This function only works if the matrix has "
522  "a single block"));
523 
524  // do a diagonal preconditioning. uses only
525  // the diagonal blocks of the matrix
526  this->block(0,0).precondition_Jacobi (dst, src, omega);
527 }
528 
529 
530 DEAL_II_NAMESPACE_CLOSE
531 
532 #endif // dealii_block_sparse_matrix_h
bool empty() const
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
virtual ~BlockSparseMatrix()
void Tvmult(BlockVector< block_number > &dst, const BlockVector< block_number > &src) const
BaseClass::value_type value_type
BlockSparseMatrix & operator=(const BlockSparseMatrix &)
const BlockSparsityPattern & get_sparsity_pattern() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1142
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:323
std::size_t memory_consumption() const
size_type n_actually_nonzero_elements(const double threshold=0.0) const
BlockMatrixBase< SparseMatrix< number > > BaseClass
static ::ExceptionBase & ExcNotQuadratic()
SmartPointer< const BlockSparsityPattern, BlockSparseMatrix< number > > sparsity_pattern
BlockSparseMatrix()=default
size_type n_nonzero_elements() const
BaseClass::BlockType BlockType
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