Reference documentation for deal.II version GIT 6a72d26406 2023-06-07 13:05:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_block_sparse_matrix.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2022 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 
17 
18 #ifdef DEAL_II_WITH_TRILINOS
19 
22 
24 
25 namespace TrilinosWrappers
26 {
28  {
29  // delete previous content of
30  // the subobjects array
31  try
32  {
33  clear();
34  }
35  catch (...)
36  {}
37  }
38 
39 
40 
41 # ifndef DOXYGEN
42  void
43  BlockSparseMatrix::reinit(const size_type n_block_rows,
44  const size_type n_block_columns)
45  {
46  // first delete previous content of
47  // the subobjects array
48  clear();
49 
50  // then resize. set sizes of blocks to
51  // zero. user will later have to call
52  // collect_sizes for this
53  this->sub_objects.reinit(n_block_rows, n_block_columns);
54  this->row_block_indices.reinit(n_block_rows, 0);
55  this->column_block_indices.reinit(n_block_columns, 0);
56 
57  // and reinitialize the blocks
58  for (size_type r = 0; r < this->n_block_rows(); ++r)
59  for (size_type c = 0; c < this->n_block_cols(); ++c)
60  {
61  BlockType *p = new BlockType();
62 
63  Assert(this->sub_objects[r][c] == nullptr, ExcInternalError());
64  this->sub_objects[r][c] = p;
65  }
66  }
67 # endif
68 
69 
70 
71  template <typename BlockSparsityPatternType>
72  void
74  const std::vector<IndexSet> & parallel_partitioning,
75  const BlockSparsityPatternType &block_sparsity_pattern,
76  const MPI_Comm communicator,
77  const bool exchange_data)
78  {
79  std::vector<Epetra_Map> epetra_maps;
80  for (size_type i = 0; i < block_sparsity_pattern.n_block_rows(); ++i)
81  epetra_maps.push_back(
82  parallel_partitioning[i].make_trilinos_map(communicator, false));
83 
84  Assert(epetra_maps.size() == block_sparsity_pattern.n_block_rows(),
85  ExcDimensionMismatch(epetra_maps.size(),
86  block_sparsity_pattern.n_block_rows()));
87  Assert(epetra_maps.size() == block_sparsity_pattern.n_block_cols(),
88  ExcDimensionMismatch(epetra_maps.size(),
89  block_sparsity_pattern.n_block_cols()));
90 
91  const size_type n_block_rows = epetra_maps.size();
92  (void)n_block_rows;
93 
94  Assert(n_block_rows == block_sparsity_pattern.n_block_rows(),
96  block_sparsity_pattern.n_block_rows()));
97  Assert(n_block_rows == block_sparsity_pattern.n_block_cols(),
99  block_sparsity_pattern.n_block_cols()));
100 
101 
102  // Call the other basic reinit function, ...
103  reinit(block_sparsity_pattern.n_block_rows(),
104  block_sparsity_pattern.n_block_cols());
105 
106  // ... set the correct sizes, ...
107  this->row_block_indices = block_sparsity_pattern.get_row_indices();
108  this->column_block_indices = block_sparsity_pattern.get_column_indices();
109 
110  // ... and then assign the correct
111  // data to the blocks.
112  for (size_type r = 0; r < this->n_block_rows(); ++r)
113  for (size_type c = 0; c < this->n_block_cols(); ++c)
114  {
115  this->sub_objects[r][c]->reinit(parallel_partitioning[r],
116  parallel_partitioning[c],
117  block_sparsity_pattern.block(r, c),
118  communicator,
119  exchange_data);
120  }
121  }
122 
123 
124 
125  template <typename BlockSparsityPatternType>
126  void
128  const BlockSparsityPatternType &block_sparsity_pattern)
129  {
130  std::vector<IndexSet> parallel_partitioning;
131  for (size_type i = 0; i < block_sparsity_pattern.n_block_rows(); ++i)
132  parallel_partitioning.emplace_back(
133  complete_index_set(block_sparsity_pattern.block(i, 0).n_rows()));
134 
135  reinit(parallel_partitioning, block_sparsity_pattern);
136  }
137 
138 
139 
140  template <>
141  void
142  BlockSparseMatrix::reinit(const BlockSparsityPattern &block_sparsity_pattern)
143  {
144  // Call the other basic reinit function, ...
145  reinit(block_sparsity_pattern.n_block_rows(),
146  block_sparsity_pattern.n_block_cols());
147 
148  // ... set the correct sizes, ...
149  this->row_block_indices = block_sparsity_pattern.get_row_indices();
150  this->column_block_indices = block_sparsity_pattern.get_column_indices();
151 
152  // ... and then assign the correct
153  // data to the blocks.
154  for (size_type r = 0; r < this->n_block_rows(); ++r)
155  for (size_type c = 0; c < this->n_block_cols(); ++c)
156  {
157  this->sub_objects[r][c]->reinit(block_sparsity_pattern.block(r, c));
158  }
159  }
160 
161 
162 
163  void
165  const std::vector<IndexSet> & parallel_partitioning,
166  const ::BlockSparseMatrix<double> &dealii_block_sparse_matrix,
167  const MPI_Comm communicator,
168  const double drop_tolerance)
169  {
170  const size_type n_block_rows = parallel_partitioning.size();
171 
172  Assert(n_block_rows == dealii_block_sparse_matrix.n_block_rows(),
174  dealii_block_sparse_matrix.n_block_rows()));
175  Assert(n_block_rows == dealii_block_sparse_matrix.n_block_cols(),
177  dealii_block_sparse_matrix.n_block_cols()));
178 
179  // Call the other basic reinit function ...
181 
182  // ... and then assign the correct
183  // data to the blocks.
184  for (size_type r = 0; r < this->n_block_rows(); ++r)
185  for (size_type c = 0; c < this->n_block_cols(); ++c)
186  {
187  this->sub_objects[r][c]->reinit(parallel_partitioning[r],
188  parallel_partitioning[c],
189  dealii_block_sparse_matrix.block(r,
190  c),
191  communicator,
192  drop_tolerance);
193  }
194 
195  collect_sizes();
196  }
197 
198 
199 
200  void
202  const ::BlockSparseMatrix<double> &dealii_block_sparse_matrix,
203  const double drop_tolerance)
204  {
205  Assert(dealii_block_sparse_matrix.n_block_rows() ==
206  dealii_block_sparse_matrix.n_block_cols(),
207  ExcDimensionMismatch(dealii_block_sparse_matrix.n_block_rows(),
208  dealii_block_sparse_matrix.n_block_cols()));
209  Assert(dealii_block_sparse_matrix.m() == dealii_block_sparse_matrix.n(),
210  ExcDimensionMismatch(dealii_block_sparse_matrix.m(),
211  dealii_block_sparse_matrix.n()));
212 
213  std::vector<IndexSet> parallel_partitioning;
214  for (size_type i = 0; i < dealii_block_sparse_matrix.n_block_rows(); ++i)
215  parallel_partitioning.emplace_back(
216  complete_index_set(dealii_block_sparse_matrix.block(i, 0).m()));
217 
218  reinit(parallel_partitioning,
219  dealii_block_sparse_matrix,
220  MPI_COMM_SELF,
221  drop_tolerance);
222  }
223 
224 
225 
226  void
228  {
229  // simply forward to the (non-public) function of the base class
231  }
232 
233 
234 
235  std::uint64_t
237  {
238  std::uint64_t n_nonzero = 0;
239  for (size_type rows = 0; rows < this->n_block_rows(); ++rows)
240  for (size_type cols = 0; cols < this->n_block_cols(); ++cols)
241  n_nonzero += this->block(rows, cols).n_nonzero_elements();
242 
243  return n_nonzero;
244  }
245 
246 
247 
250  const MPI::BlockVector &x,
251  const MPI::BlockVector &b) const
252  {
253  vmult(dst, x);
254  dst -= b;
255  dst *= -1.;
256 
257  return dst.l2_norm();
258  }
259 
260 
261 
262  // TODO: In the following we
263  // use the same code as just
264  // above three more times. Use
265  // templates.
268  const MPI::Vector & x,
269  const MPI::BlockVector &b) const
270  {
271  vmult(dst, x);
272  dst -= b;
273  dst *= -1.;
274 
275  return dst.l2_norm();
276  }
277 
278 
279 
282  const MPI::BlockVector &x,
283  const MPI::Vector & b) const
284  {
285  vmult(dst, x);
286  dst -= b;
287  dst *= -1.;
288 
289  return dst.l2_norm();
290  }
291 
292 
293 
296  const MPI::Vector &x,
297  const MPI::Vector &b) const
298  {
299  vmult(dst, x);
300  dst -= b;
301  dst *= -1.;
302 
303  return dst.l2_norm();
304  }
305 
306 
307 
308  MPI_Comm
310  {
311  Assert(this->n_block_cols() != 0, ExcNotInitialized());
312  Assert(this->n_block_rows() != 0, ExcNotInitialized());
313  return this->sub_objects[0][0]->get_mpi_communicator();
314  }
315 
316 
317 
318 # ifndef DOXYGEN
319  // -------------------- explicit instantiations -----------------------
320  //
321  template void
322  BlockSparseMatrix::reinit(const ::BlockSparsityPattern &);
323  template void
324  BlockSparseMatrix::reinit(const ::BlockDynamicSparsityPattern &);
325 
326  template void
327  BlockSparseMatrix::reinit(const std::vector<IndexSet> &,
328  const ::BlockDynamicSparsityPattern &,
329  const MPI_Comm,
330  const bool);
331 # endif // DOXYGEN
332 
333 } // namespace TrilinosWrappers
334 
335 
337 
338 #endif
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
BlockType & block(const unsigned int row, const unsigned int column)
unsigned int n_block_rows() const
Table< 2, SmartPointer< BlockType, BlockMatrixBase< SparseMatrix > > > sub_objects
unsigned int n_block_cols() const
SparsityPatternType & block(const size_type row, const size_type column)
const BlockIndices & get_column_indices() const
const BlockIndices & get_row_indices() const
real_type l2_norm() const
std::size_t n_nonzero_elements() const
TrilinosScalar residual(MPI::BlockVector &dst, const MPI::BlockVector &x, const MPI::BlockVector &b) const
void reinit(const size_type n_block_rows, const size_type n_block_columns)
void vmult(VectorType1 &dst, const VectorType2 &src) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1614
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1089
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)