Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\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
37 # ifndef DOXYGEN
39 template <typename number>
40 class BlockSparseMatrix;
41 # endif
42 
43 namespace TrilinosWrappers
44 {
72  class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix>
73  {
74  public:
79 
84 
89  using pointer = BaseClass::pointer;
90  using const_pointer = BaseClass::const_pointer;
91  using reference = BaseClass::reference;
92  using const_reference = BaseClass::const_reference;
93  using size_type = BaseClass::size_type;
96 
108  BlockSparseMatrix() = default;
109 
113  ~BlockSparseMatrix() override;
114 
120  operator=(const BlockSparseMatrix &) = default;
121 
132  operator=(const double d);
133 
147  void
148  reinit(const size_type n_block_rows, const size_type n_block_columns);
149 
155  template <typename BlockSparsityPatternType>
156  void
157  reinit(const std::vector<Epetra_Map> & input_maps,
158  const BlockSparsityPatternType &block_sparsity_pattern,
159  const bool exchange_data = false);
160 
166  template <typename BlockSparsityPatternType>
167  void
168  reinit(const std::vector<IndexSet> & input_maps,
169  const BlockSparsityPatternType &block_sparsity_pattern,
170  const MPI_Comm & communicator = MPI_COMM_WORLD,
171  const bool exchange_data = false);
172 
178  template <typename BlockSparsityPatternType>
179  void
180  reinit(const BlockSparsityPatternType &block_sparsity_pattern);
181 
190  DEAL_II_DEPRECATED
191  void
192  reinit(const std::vector<Epetra_Map> & input_maps,
193  const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
194  const double drop_tolerance = 1e-13);
195 
203  void
204  reinit(const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
205  const double drop_tolerance = 1e-13);
206 
214  bool
215  is_compressed() const;
216 
227  void
228  collect_sizes();
229 
234  size_type
235  n_nonzero_elements() const;
236 
240  MPI_Comm
241  get_mpi_communicator() const;
242 
252  DEAL_II_DEPRECATED
253  std::vector<Epetra_Map>
254  domain_partitioner() const;
255 
265  DEAL_II_DEPRECATED
266  std::vector<Epetra_Map>
267  range_partitioner() const;
268 
274  std::vector<IndexSet>
276 
282  std::vector<IndexSet>
284 
291  template <typename VectorType1, typename VectorType2>
292  void
293  vmult(VectorType1 &dst, const VectorType2 &src) const;
294 
300  template <typename VectorType1, typename VectorType2>
301  void
302  Tvmult(VectorType1 &dst, const VectorType2 &src) const;
303 
316  TrilinosScalar
318  const MPI::BlockVector &x,
319  const MPI::BlockVector &b) const;
320 
328  TrilinosScalar
330  const MPI::Vector & x,
331  const MPI::BlockVector &b) const;
332 
340  TrilinosScalar
341  residual(MPI::Vector & dst,
342  const MPI::BlockVector &x,
343  const MPI::Vector & b) const;
344 
352  TrilinosScalar
353  residual(MPI::Vector & dst,
354  const MPI::Vector &x,
355  const MPI::Vector &b) const;
356 
362 
372  int,
373  int,
374  int,
375  int,
376  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
377  << ',' << arg4 << "] have differing row numbers.");
378 
383  int,
384  int,
385  int,
386  int,
387  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
388  << ',' << arg4 << "] have differing column numbers.");
390 
391  private:
395  template <typename VectorType1, typename VectorType2>
396  void
397  vmult(VectorType1 & dst,
398  const VectorType2 &src,
399  const bool transpose,
400  const std::integral_constant<bool, true>,
401  const std::integral_constant<bool, true>) const;
402 
407  template <typename VectorType1, typename VectorType2>
408  void
409  vmult(VectorType1 & dst,
410  const VectorType2 &src,
411  const bool transpose,
412  const std::integral_constant<bool, false>,
413  const std::integral_constant<bool, true>) const;
414 
419  template <typename VectorType1, typename VectorType2>
420  void
421  vmult(VectorType1 & dst,
422  const VectorType2 &src,
423  const bool transpose,
424  const std::integral_constant<bool, true>,
425  const std::integral_constant<bool, false>) const;
426 
432  template <typename VectorType1, typename VectorType2>
433  void
434  vmult(VectorType1 & dst,
435  const VectorType2 &src,
436  const bool transpose,
437  const std::integral_constant<bool, false>,
438  const std::integral_constant<bool, false>) const;
439  };
440 
441 
442 
445  // ------------- inline and template functions -----------------
446 
447 
448 
449  inline BlockSparseMatrix &
451  {
453 
454  for (size_type r = 0; r < this->n_block_rows(); ++r)
455  for (size_type c = 0; c < this->n_block_cols(); ++c)
456  this->block(r, c) = d;
457 
458  return *this;
459  }
460 
461 
462 
463  inline bool
465  {
466  bool compressed = true;
467  for (size_type row = 0; row < n_block_rows(); ++row)
468  for (size_type col = 0; col < n_block_cols(); ++col)
469  if (block(row, col).is_compressed() == false)
470  {
471  compressed = false;
472  break;
473  }
474 
475  return compressed;
476  }
477 
478 
479 
480  template <typename VectorType1, typename VectorType2>
481  inline void
482  BlockSparseMatrix::vmult(VectorType1 &dst, const VectorType2 &src) const
483  {
484  vmult(dst,
485  src,
486  false,
487  std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
488  std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
489  }
490 
491 
492 
493  template <typename VectorType1, typename VectorType2>
494  inline void
495  BlockSparseMatrix::Tvmult(VectorType1 &dst, const VectorType2 &src) const
496  {
497  vmult(dst,
498  src,
499  true,
500  std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
501  std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
502  }
503 
504 
505 
506  template <typename VectorType1, typename VectorType2>
507  inline void
508  BlockSparseMatrix::vmult(VectorType1 & dst,
509  const VectorType2 &src,
510  const bool transpose,
511  std::integral_constant<bool, true>,
512  std::integral_constant<bool, true>) const
513  {
514  if (transpose == true)
516  else
518  }
519 
520 
521 
522  template <typename VectorType1, typename VectorType2>
523  inline void
524  BlockSparseMatrix::vmult(VectorType1 & dst,
525  const VectorType2 &src,
526  const bool transpose,
527  std::integral_constant<bool, false>,
528  std::integral_constant<bool, true>) const
529  {
530  if (transpose == true)
532  else
534  }
535 
536 
537 
538  template <typename VectorType1, typename VectorType2>
539  inline void
540  BlockSparseMatrix::vmult(VectorType1 & dst,
541  const VectorType2 &src,
542  const bool transpose,
543  std::integral_constant<bool, true>,
544  std::integral_constant<bool, false>) const
545  {
546  if (transpose == true)
548  else
550  }
551 
552 
553 
554  template <typename VectorType1, typename VectorType2>
555  inline void
556  BlockSparseMatrix::vmult(VectorType1 & dst,
557  const VectorType2 &src,
558  const bool transpose,
559  std::integral_constant<bool, false>,
560  std::integral_constant<bool, false>) const
561  {
562  if (transpose == true)
564  else
566  }
567 
568 
569 
570  inline std::vector<IndexSet>
572  {
573  Assert(this->n_block_cols() != 0, ExcNotInitialized());
574  Assert(this->n_block_rows() != 0, ExcNotInitialized());
575 
576  std::vector<IndexSet> domain_indices;
577  for (size_type c = 0; c < this->n_block_cols(); ++c)
578  domain_indices.push_back(
579  this->sub_objects[0][c]->locally_owned_domain_indices());
580 
581  return domain_indices;
582  }
583 
584 
585 
586  inline std::vector<IndexSet>
588  {
589  Assert(this->n_block_cols() != 0, ExcNotInitialized());
590  Assert(this->n_block_rows() != 0, ExcNotInitialized());
591 
592  std::vector<IndexSet> range_indices;
593  for (size_type r = 0; r < this->n_block_rows(); ++r)
594  range_indices.push_back(
595  this->sub_objects[r][0]->locally_owned_range_indices());
596 
597  return range_indices;
598  }
599 
600 
601 
602  namespace internal
603  {
604  namespace BlockLinearOperatorImplementation
605  {
621  template <typename PayloadBlockType>
623  {
624  public:
628  using BlockType = PayloadBlockType;
629 
639  template <typename... Args>
640  TrilinosBlockPayload(const Args &...)
641  {
642  static_assert(
643  std::is_same<
644  PayloadBlockType,
646  "TrilinosBlockPayload can only accept a payload of type TrilinosPayload.");
647  }
648  };
649 
650  } // namespace BlockLinearOperatorImplementation
651  } /* namespace internal */
652 
653 
654 } /* namespace TrilinosWrappers */
655 
656 
657 DEAL_II_NAMESPACE_CLOSE
658 
659 #endif // DEAL_II_WITH_TRILINOS
660 
661 #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
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
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)
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