1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2019 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_mg_matrix_h
17 #define dealii_mg_matrix_h
18
19 #include <deal.II/base/config.h>
20
21 #include <deal.II/base/mg_level_object.h>
22
23 #include <deal.II/lac/linear_operator.h>
24 #include <deal.II/lac/sparse_matrix.h>
25 #include <deal.II/lac/vector.h>
26
27 #include <deal.II/multigrid/mg_base.h>
28
29 #include <memory>
30
31 DEAL_II_NAMESPACE_OPEN
32
33 /*!@addtogroup mg */
34 /*@{*/
35
36 namespace mg
37 {
38 /**
39 * Multilevel matrix. This matrix stores an MGLevelObject of
40 * LinearOperator objects. It implements the interface defined in
41 * MGMatrixBase, so that it can be used as a matrix in Multigrid.
42 */
43 template <typename VectorType = Vector<double>>
44 class Matrix : public MGMatrixBase<VectorType>
45 {
46 public:
47 /**
48 * Default constructor for an empty object.
49 */
50 Matrix() = default;
51
52 /**
53 * Constructor setting up pointers to the matrices in <tt>M</tt> by
54 * calling initialize().
55 */
56 template <typename MatrixType>
57 Matrix(const MGLevelObject<MatrixType> &M);
58
59 /**
60 * Initialize the object such that the level multiplication uses the
61 * matrices in <tt>M</tt>
62 */
63 template <typename MatrixType>
64 void
65 initialize(const MGLevelObject<MatrixType> &M);
66
67 /**
68 * Reset the object.
69 */
70 void
71 reset();
72
73 /**
74 * Access matrix on a level.
75 */
76 const LinearOperator<VectorType> &operator[](unsigned int level) const;
77
78 virtual void
79 vmult(const unsigned int level,
80 VectorType & dst,
81 const VectorType & src) const override;
82 virtual void
83 vmult_add(const unsigned int level,
84 VectorType & dst,
85 const VectorType & src) const override;
86 virtual void
87 Tvmult(const unsigned int level,
88 VectorType & dst,
89 const VectorType & src) const override;
90 virtual void
91 Tvmult_add(const unsigned int level,
92 VectorType & dst,
93 const VectorType & src) const override;
94 virtual unsigned int
95 get_minlevel() const override;
96 virtual unsigned int
97 get_maxlevel() const override;
98
99 /**
100 * Memory used by this object.
101 */
102 std::size_t
103 memory_consumption() const;
104
105 private:
106 MGLevelObject<LinearOperator<VectorType>> matrices;
107 };
108
109 } // namespace mg
110
111
112 /**
113 * Multilevel matrix selecting from block matrices. This class implements the
114 * interface defined by MGMatrixBase. The template parameter @p MatrixType
115 * should be a block matrix class like BlockSparseMatrix or @p
116 * BlockSparseMatrixEZ. Then, this class stores a pointer to a MGLevelObject
117 * of this matrix class. In each @p vmult, the block selected on
118 * initialization will be multiplied with the vector provided.
119 */
120 template <typename MatrixType, typename number>
121 class MGMatrixSelect : public MGMatrixBase<Vector<number>>
122 {
123 public:
124 /**
125 * Constructor. @p row and @p col are the coordinate of the selected block.
126 * The other argument is handed over to the @p SmartPointer constructor.
127 */
128 MGMatrixSelect(const unsigned int row = 0,
129 const unsigned int col = 0,
130 MGLevelObject<MatrixType> *matrix = 0);
131
132 /**
133 * Set the matrix object to be used. The matrix object must exist longer as
134 * the @p MGMatrixSelect object, since only a pointer is stored.
135 */
136 void
137 set_matrix(MGLevelObject<MatrixType> *M);
138
139 /**
140 * Select the block for multiplication.
141 */
142 void
143 select_block(const unsigned int row, const unsigned int col);
144
145 /**
146 * Matrix-vector-multiplication on a certain level.
147 */
148 virtual void
149 vmult(const unsigned int level,
150 Vector<number> & dst,
151 const Vector<number> &src) const;
152
153 /**
154 * Adding matrix-vector-multiplication on a certain level.
155 */
156 virtual void
157 vmult_add(const unsigned int level,
158 Vector<number> & dst,
159 const Vector<number> &src) const;
160
161 /**
162 * Transpose matrix-vector-multiplication on a certain level.
163 */
164 virtual void
165 Tvmult(const unsigned int level,
166 Vector<number> & dst,
167 const Vector<number> &src) const;
168
169 /**
170 * Adding transpose matrix-vector-multiplication on a certain level.
171 */
172 virtual void
173 Tvmult_add(const unsigned int level,
174 Vector<number> & dst,
175 const Vector<number> &src) const;
176
177 private:
178 /**
179 * Pointer to the matrix objects on each level.
180 */
181 SmartPointer<MGLevelObject<MatrixType>, MGMatrixSelect<MatrixType, number>>
182 matrix;
183 /**
184 * Row coordinate of selected block.
185 */
186 unsigned int row;
187 /**
188 * Column coordinate of selected block.
189 */
190 unsigned int col;
191 };
192
193 /*@}*/
194
195 /*----------------------------------------------------------------------*/
196
197 namespace mg
198 {
199 template <typename VectorType>
200 template <typename MatrixType>
201 inline void
initialize(const MGLevelObject<MatrixType> & p)202 Matrix<VectorType>::initialize(const MGLevelObject<MatrixType> &p)
203 {
204 matrices.resize(p.min_level(), p.max_level());
205 for (unsigned int level = p.min_level(); level <= p.max_level(); ++level)
206 {
207 // Workaround: Unfortunately, not every "p[level]" object has a
208 // rich enough interface to populate reinit_(domain|range)_vector.
209 // Thus, apply an empty LinearOperator exemplar.
210 matrices[level] =
211 linear_operator<VectorType>(LinearOperator<VectorType>(), p[level]);
212 }
213 }
214
215
216
217 template <typename VectorType>
218 inline void
reset()219 Matrix<VectorType>::reset()
220 {
221 matrices.resize(0, 0);
222 }
223
224
225
226 template <typename VectorType>
227 template <typename MatrixType>
Matrix(const MGLevelObject<MatrixType> & p)228 inline Matrix<VectorType>::Matrix(const MGLevelObject<MatrixType> &p)
229 {
230 initialize(p);
231 }
232
233
234
235 template <typename VectorType>
236 inline const LinearOperator<VectorType> &Matrix<VectorType>::
237 operator[](unsigned int level) const
238 {
239 return matrices[level];
240 }
241
242
243
244 template <typename VectorType>
245 void
vmult(const unsigned int level,VectorType & dst,const VectorType & src)246 Matrix<VectorType>::vmult(const unsigned int level,
247 VectorType & dst,
248 const VectorType & src) const
249 {
250 matrices[level].vmult(dst, src);
251 }
252
253
254
255 template <typename VectorType>
256 void
vmult_add(const unsigned int level,VectorType & dst,const VectorType & src)257 Matrix<VectorType>::vmult_add(const unsigned int level,
258 VectorType & dst,
259 const VectorType & src) const
260 {
261 matrices[level].vmult_add(dst, src);
262 }
263
264
265
266 template <typename VectorType>
267 void
Tvmult(const unsigned int level,VectorType & dst,const VectorType & src)268 Matrix<VectorType>::Tvmult(const unsigned int level,
269 VectorType & dst,
270 const VectorType & src) const
271 {
272 matrices[level].Tvmult(dst, src);
273 }
274
275
276
277 template <typename VectorType>
278 void
Tvmult_add(const unsigned int level,VectorType & dst,const VectorType & src)279 Matrix<VectorType>::Tvmult_add(const unsigned int level,
280 VectorType & dst,
281 const VectorType & src) const
282 {
283 matrices[level].Tvmult_add(dst, src);
284 }
285
286
287
288 template <typename VectorType>
289 unsigned int
get_minlevel()290 Matrix<VectorType>::get_minlevel() const
291 {
292 return matrices.min_level();
293 }
294
295
296
297 template <typename VectorType>
298 unsigned int
get_maxlevel()299 Matrix<VectorType>::get_maxlevel() const
300 {
301 return matrices.max_level();
302 }
303
304
305
306 template <typename VectorType>
307 inline std::size_t
memory_consumption()308 Matrix<VectorType>::memory_consumption() const
309 {
310 return sizeof(*this) + matrices->memory_consumption();
311 }
312 } // namespace mg
313
314
315 /*----------------------------------------------------------------------*/
316
317 template <typename MatrixType, typename number>
MGMatrixSelect(const unsigned int row,const unsigned int col,MGLevelObject<MatrixType> * p)318 MGMatrixSelect<MatrixType, number>::MGMatrixSelect(const unsigned int row,
319 const unsigned int col,
320 MGLevelObject<MatrixType> *p)
321 : matrix(p, typeid(*this).name())
322 , row(row)
323 , col(col)
324 {}
325
326
327
328 template <typename MatrixType, typename number>
329 void
set_matrix(MGLevelObject<MatrixType> * p)330 MGMatrixSelect<MatrixType, number>::set_matrix(MGLevelObject<MatrixType> *p)
331 {
332 matrix = p;
333 }
334
335
336
337 template <typename MatrixType, typename number>
338 void
select_block(const unsigned int brow,const unsigned int bcol)339 MGMatrixSelect<MatrixType, number>::select_block(const unsigned int brow,
340 const unsigned int bcol)
341 {
342 row = brow;
343 col = bcol;
344 }
345
346
347
348 template <typename MatrixType, typename number>
349 void
vmult(const unsigned int level,Vector<number> & dst,const Vector<number> & src)350 MGMatrixSelect<MatrixType, number>::vmult(const unsigned int level,
351 Vector<number> & dst,
352 const Vector<number> &src) const
353 {
354 Assert(matrix != 0, ExcNotInitialized());
355
356 const MGLevelObject<MatrixType> &m = *matrix;
357 m[level].block(row, col).vmult(dst, src);
358 }
359
360
361
362 template <typename MatrixType, typename number>
363 void
vmult_add(const unsigned int level,Vector<number> & dst,const Vector<number> & src)364 MGMatrixSelect<MatrixType, number>::vmult_add(const unsigned int level,
365 Vector<number> & dst,
366 const Vector<number> &src) const
367 {
368 Assert(matrix != 0, ExcNotInitialized());
369
370 const MGLevelObject<MatrixType> &m = *matrix;
371 m[level].block(row, col).vmult_add(dst, src);
372 }
373
374
375
376 template <typename MatrixType, typename number>
377 void
Tvmult(const unsigned int level,Vector<number> & dst,const Vector<number> & src)378 MGMatrixSelect<MatrixType, number>::Tvmult(const unsigned int level,
379 Vector<number> & dst,
380 const Vector<number> &src) const
381 {
382 Assert(matrix != 0, ExcNotInitialized());
383
384 const MGLevelObject<MatrixType> &m = *matrix;
385 m[level].block(row, col).Tvmult(dst, src);
386 }
387
388
389
390 template <typename MatrixType, typename number>
391 void
Tvmult_add(const unsigned int level,Vector<number> & dst,const Vector<number> & src)392 MGMatrixSelect<MatrixType, number>::Tvmult_add(const unsigned int level,
393 Vector<number> & dst,
394 const Vector<number> &src) const
395 {
396 Assert(matrix != 0, ExcNotInitialized());
397
398 const MGLevelObject<MatrixType> &m = *matrix;
399 m[level].block(row, col).Tvmult_add(dst, src);
400 }
401
402 DEAL_II_NAMESPACE_CLOSE
403
404 #endif
405