Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
trilinos_block_sparse_matrix.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2023 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
25namespace 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
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
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)
unsigned int n_block_rows() const
Table< 2, SmartPointer< BlockType, BlockMatrixBase< SparseMatrix > > > sub_objects
unsigned int n_block_cols() const
BlockType & block(const unsigned int row, const unsigned int column)
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1089
double TrilinosScalar
Definition types.h:175