Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
trilinos_block_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 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_trilinos_block_sparse_matrix_h
17 #define dealii_trilinos_block_sparse_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/template_constraints.h>
25 
26 # include <deal.II/lac/block_matrix_base.h>
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/full_matrix.h>
29 # include <deal.II/lac/trilinos_parallel_block_vector.h>
30 # include <deal.II/lac/trilinos_sparse_matrix.h>
31 
32 # include <cmath>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
36 // forward declarations
38 template <typename number>
39 class BlockSparseMatrix;
40 
41 
42 namespace TrilinosWrappers
43 {
71  class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix>
72  {
73  public:
78 
83 
88  using pointer = BaseClass::pointer;
89  using const_pointer = BaseClass::const_pointer;
90  using reference = BaseClass::reference;
91  using const_reference = BaseClass::const_reference;
92  using size_type = BaseClass::size_type;
95 
107  BlockSparseMatrix() = default;
108 
112  ~BlockSparseMatrix() override;
113 
119  operator=(const BlockSparseMatrix &) = default;
120 
131  operator=(const double d);
132 
146  void
147  reinit(const size_type n_block_rows, const size_type n_block_columns);
148 
154  template <typename BlockSparsityPatternType>
155  void
156  reinit(const std::vector<Epetra_Map> & input_maps,
157  const BlockSparsityPatternType &block_sparsity_pattern,
158  const bool exchange_data = false);
159 
165  template <typename BlockSparsityPatternType>
166  void
167  reinit(const std::vector<IndexSet> & input_maps,
168  const BlockSparsityPatternType &block_sparsity_pattern,
169  const MPI_Comm & communicator = MPI_COMM_WORLD,
170  const bool exchange_data = false);
171 
177  template <typename BlockSparsityPatternType>
178  void
179  reinit(const BlockSparsityPatternType &block_sparsity_pattern);
180 
189  DEAL_II_DEPRECATED
190  void
191  reinit(const std::vector<Epetra_Map> & input_maps,
192  const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
193  const double drop_tolerance = 1e-13);
194 
202  void
203  reinit(const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
204  const double drop_tolerance = 1e-13);
205 
213  bool
214  is_compressed() const;
215 
226  void
227  collect_sizes();
228 
233  size_type
234  n_nonzero_elements() const;
235 
239  MPI_Comm
240  get_mpi_communicator() const;
241 
251  DEAL_II_DEPRECATED
252  std::vector<Epetra_Map>
253  domain_partitioner() const;
254 
264  DEAL_II_DEPRECATED
265  std::vector<Epetra_Map>
266  range_partitioner() const;
267 
273  std::vector<IndexSet>
275 
281  std::vector<IndexSet>
283 
290  template <typename VectorType1, typename VectorType2>
291  void
292  vmult(VectorType1 &dst, const VectorType2 &src) const;
293 
299  template <typename VectorType1, typename VectorType2>
300  void
301  Tvmult(VectorType1 &dst, const VectorType2 &src) const;
302 
315  TrilinosScalar
317  const MPI::BlockVector &x,
318  const MPI::BlockVector &b) const;
319 
327  TrilinosScalar
329  const MPI::Vector & x,
330  const MPI::BlockVector &b) const;
331 
339  TrilinosScalar
340  residual(MPI::Vector & dst,
341  const MPI::BlockVector &x,
342  const MPI::Vector & b) const;
343 
351  TrilinosScalar
352  residual(MPI::Vector & dst,
353  const MPI::Vector &x,
354  const MPI::Vector &b) const;
355 
361 
371  int,
372  int,
373  int,
374  int,
375  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
376  << ',' << arg4 << "] have differing row numbers.");
377 
382  int,
383  int,
384  int,
385  int,
386  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
387  << ',' << arg4 << "] have differing column numbers.");
389 
390  private:
394  template <typename VectorType1, typename VectorType2>
395  void
396  vmult(VectorType1 & dst,
397  const VectorType2 &src,
398  const bool transpose,
399  const std::integral_constant<bool, true>,
400  const std::integral_constant<bool, true>) const;
401 
406  template <typename VectorType1, typename VectorType2>
407  void
408  vmult(VectorType1 & dst,
409  const VectorType2 &src,
410  const bool transpose,
411  const std::integral_constant<bool, false>,
412  const std::integral_constant<bool, true>) const;
413 
418  template <typename VectorType1, typename VectorType2>
419  void
420  vmult(VectorType1 & dst,
421  const VectorType2 &src,
422  const bool transpose,
423  const std::integral_constant<bool, true>,
424  const std::integral_constant<bool, false>) const;
425 
431  template <typename VectorType1, typename VectorType2>
432  void
433  vmult(VectorType1 & dst,
434  const VectorType2 &src,
435  const bool transpose,
436  const std::integral_constant<bool, false>,
437  const std::integral_constant<bool, false>) const;
438  };
439 
440 
441 
444  // ------------- inline and template functions -----------------
445 
446 
447 
448  inline BlockSparseMatrix &
450  {
452 
453  for (size_type r = 0; r < this->n_block_rows(); ++r)
454  for (size_type c = 0; c < this->n_block_cols(); ++c)
455  this->block(r, c) = d;
456 
457  return *this;
458  }
459 
460 
461 
462  inline bool
464  {
465  bool compressed = true;
466  for (size_type row = 0; row < n_block_rows(); ++row)
467  for (size_type col = 0; col < n_block_cols(); ++col)
468  if (block(row, col).is_compressed() == false)
469  {
470  compressed = false;
471  break;
472  }
473 
474  return compressed;
475  }
476 
477 
478 
479  template <typename VectorType1, typename VectorType2>
480  inline void
481  BlockSparseMatrix::vmult(VectorType1 &dst, const VectorType2 &src) const
482  {
483  vmult(dst,
484  src,
485  false,
486  std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
487  std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
488  }
489 
490 
491 
492  template <typename VectorType1, typename VectorType2>
493  inline void
494  BlockSparseMatrix::Tvmult(VectorType1 &dst, const VectorType2 &src) const
495  {
496  vmult(dst,
497  src,
498  true,
499  std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
500  std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
501  }
502 
503 
504 
505  template <typename VectorType1, typename VectorType2>
506  inline void
507  BlockSparseMatrix::vmult(VectorType1 & dst,
508  const VectorType2 &src,
509  const bool transpose,
510  std::integral_constant<bool, true>,
511  std::integral_constant<bool, true>) const
512  {
513  if (transpose == true)
515  else
517  }
518 
519 
520 
521  template <typename VectorType1, typename VectorType2>
522  inline void
523  BlockSparseMatrix::vmult(VectorType1 & dst,
524  const VectorType2 &src,
525  const bool transpose,
526  std::integral_constant<bool, false>,
527  std::integral_constant<bool, true>) const
528  {
529  if (transpose == true)
531  else
533  }
534 
535 
536 
537  template <typename VectorType1, typename VectorType2>
538  inline void
539  BlockSparseMatrix::vmult(VectorType1 & dst,
540  const VectorType2 &src,
541  const bool transpose,
542  std::integral_constant<bool, true>,
543  std::integral_constant<bool, false>) const
544  {
545  if (transpose == true)
547  else
549  }
550 
551 
552 
553  template <typename VectorType1, typename VectorType2>
554  inline void
555  BlockSparseMatrix::vmult(VectorType1 & dst,
556  const VectorType2 &src,
557  const bool transpose,
558  std::integral_constant<bool, false>,
559  std::integral_constant<bool, false>) const
560  {
561  if (transpose == true)
563  else
565  }
566 
567 
568 
569  inline std::vector<IndexSet>
571  {
572  Assert(this->n_block_cols() != 0, ExcNotInitialized());
573  Assert(this->n_block_rows() != 0, ExcNotInitialized());
574 
575  std::vector<IndexSet> domain_indices;
576  for (size_type c = 0; c < this->n_block_cols(); ++c)
577  domain_indices.push_back(
578  this->sub_objects[0][c]->locally_owned_domain_indices());
579 
580  return domain_indices;
581  }
582 
583 
584 
585  inline std::vector<IndexSet>
587  {
588  Assert(this->n_block_cols() != 0, ExcNotInitialized());
589  Assert(this->n_block_rows() != 0, ExcNotInitialized());
590 
591  std::vector<IndexSet> range_indices;
592  for (size_type r = 0; r < this->n_block_rows(); ++r)
593  range_indices.push_back(
594  this->sub_objects[r][0]->locally_owned_range_indices());
595 
596  return range_indices;
597  }
598 
599 
600 
601  namespace internal
602  {
603  namespace BlockLinearOperatorImplementation
604  {
620  template <typename PayloadBlockType>
622  {
623  public:
627  using BlockType = PayloadBlockType;
628 
638  template <typename... Args>
639  TrilinosBlockPayload(const Args &...)
640  {
641  static_assert(
642  std::is_same<
643  PayloadBlockType,
645  "TrilinosBlockPayload can only accept a payload of type TrilinosPayload.");
646  }
647  };
648 
649  } // namespace BlockLinearOperatorImplementation
650  } /* namespace internal */
651 
652 
653 } /* namespace TrilinosWrappers */
654 
655 
656 DEAL_II_NAMESPACE_CLOSE
657 
658 #endif // DEAL_II_WITH_TRILINOS
659 
660 #endif // dealii_trilinos_block_sparse_matrix_h
void vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
std::vector< Epetra_Map > range_partitioner() const
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
static ::ExceptionBase & ExcNotInitialized()
void Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const
void Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
void Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const
void vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const
typename BlockType::value_type value_type
void Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
unsigned int n_block_cols() const
std::vector< Epetra_Map > domain_partitioner() const
void vmult(VectorType1 &dst, const VectorType2 &src) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
void reinit(const size_type n_block_rows, const size_type n_block_columns)
SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
void Tvmult(VectorType1 &dst, const VectorType2 &src) const
void vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const
BlockType & block(const unsigned int row, const unsigned int column)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
unsigned int n_block_rows() const
std::vector< IndexSet > locally_owned_range_indices() const
TrilinosScalar residual(MPI::BlockVector &dst, const MPI::BlockVector &x, const MPI::BlockVector &b) const
BlockSparseMatrix & operator=(const BlockSparseMatrix &)=default
std::vector< IndexSet > locally_owned_domain_indices() const