MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
blockoperator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_BLOCKOPERATOR
13 #define MFEM_BLOCKOPERATOR
14 
15 #include "../config/config.hpp"
16 #include "../general/array.hpp"
17 #include "operator.hpp"
18 #include "blockvector.hpp"
19 
20 namespace mfem
21 {
22 
23 //! @class BlockOperator
24 /**
25  * \brief A class to handle Block systems in a matrix-free implementation.
26  *
27  * Usage:
28  * - Use one of the constructors to define the block structure.
29  * - Use SetDiagonalBlock or SetBlock to fill the BlockOperator
30  * - Use the method Mult and MultTranspose to apply the operator to a vector.
31  *
32  * If a block is not set, it is assumed to be a zero block.
33  */
34 class BlockOperator : public Operator
35 {
36 public:
37  //! Constructor for BlockOperators with the same block-structure for rows and
38  //! columns.
39  /**
40  * offsets: offsets that mark the start of each row/column block (size
41  * nRowBlocks+1). Note: BlockOperator will not own/copy the data contained
42  * in offsets.
43  */
44  BlockOperator(const Array<int> & offsets);
45  //! Constructor for general BlockOperators.
46  /**
47  * row_offsets: offsets that mark the start of each row block (size
48  * nRowBlocks+1). col_offsets: offsets that mark the start of each column
49  * block (size nColBlocks+1). Note: BlockOperator will not own/copy the
50  * data contained in offsets.
51  */
52  BlockOperator(const Array<int> & row_offsets, const Array<int> & col_offsets);
53 
54  /// Copy assignment is not supported
55  BlockOperator &operator=(const BlockOperator &) = delete;
56 
57  /// Move assignment is not supported
59 
60  //! Add block op in the block-entry (iblock, iblock).
61  /**
62  * iblock: The block will be inserted in location (iblock, iblock).
63  * op: the Operator to be inserted.
64  * c: optional scalar multiple for this block.
65  */
66  void SetDiagonalBlock(int iblock, Operator *op, double c = 1.0);
67  //! Add a block op in the block-entry (iblock, jblock).
68  /**
69  * irow, icol: The block will be inserted in location (irow, icol).
70  * op: the Operator to be inserted.
71  * c: optional scalar multiple for this block.
72  */
73  void SetBlock(int iRow, int iCol, Operator *op, double c = 1.0);
74 
75  //! Return the number of row blocks
76  int NumRowBlocks() const { return nRowBlocks; }
77  //! Return the number of column blocks
78  int NumColBlocks() const { return nColBlocks; }
79 
80  //! Check if block (i,j) is a zero block
81  int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; }
82  //! Return a reference to block i,j
83  Operator & GetBlock(int i, int j)
84  { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
85  //! Return a reference to block i,j (const version)
86  const Operator & GetBlock(int i, int j) const
87  { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
88  //! Return the coefficient for block i,j
89  double GetBlockCoef(int i, int j) const
90  { MFEM_VERIFY(op(i,j), ""); return coef(i,j); }
91  //! Set the coefficient for block i,j
92  void SetBlockCoef(int i, int j, double c)
93  { MFEM_VERIFY(op(i,j), ""); coef(i,j) = c; }
94 
95  //! Return the row offsets for block starts
96  Array<int> & RowOffsets() { return row_offsets; }
97  //! Read only access to the row offsets for block starts
98  const Array<int> & RowOffsets() const { return row_offsets; }
99  //! Return the columns offsets for block starts
100  Array<int> & ColOffsets() { return col_offsets; }
101  //! Read only access to the columns offsets for block starts
102  const Array<int> & ColOffsets() const { return col_offsets; }
103 
104  /// Operator application
105  virtual void Mult (const Vector & x, Vector & y) const;
106 
107  /// Action of the transpose operator
108  virtual void MultTranspose (const Vector & x, Vector & y) const;
109 
110  ~BlockOperator();
111 
112  //! Controls the ownership of the blocks: if nonzero, BlockOperator will
113  //! delete all blocks that are set (non-NULL); the default value is zero.
115 
116 private:
117  //! Number of block rows
118  int nRowBlocks;
119  //! Number of block columns
120  int nColBlocks;
121  //! Row offsets for the starting position of each block
122  Array<int> row_offsets;
123  //! Column offsets for the starting position of each block
124  Array<int> col_offsets;
125  //! 2D array that stores each block of the operator.
127  //! 2D array that stores a coefficient for each block of the operator.
128  Array2D<double> coef;
129 
130  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose methods.
131  mutable BlockVector xblock;
132  mutable BlockVector yblock;
133  mutable Vector tmp;
134 };
135 
136 //! @class BlockDiagonalPreconditioner
137 /**
138  * \brief A class to handle Block diagonal preconditioners in a matrix-free implementation.
139  *
140  * Usage:
141  * - Use the constructors to define the block structure
142  * - Use SetDiagonalBlock to fill the BlockOperator
143  * - Use the method Mult and MultTranspose to apply the operator to a vector.
144  *
145  * If a block is not set, it is assumed to be an identity block.
146  *
147  */
149 {
150 public:
151  //! Constructor that specifies the block structure
152  BlockDiagonalPreconditioner(const Array<int> & offsets);
153  //! Add a square block op in the block-entry (iblock, iblock).
154  /**
155  * iblock: The block will be inserted in location (iblock, iblock).
156  * op: the Operator to be inserted.
157  */
158  void SetDiagonalBlock(int iblock, Operator *op);
159  //! This method is present since required by the abstract base class Solver
160  virtual void SetOperator(const Operator &op) { }
161 
162  //! Return the number of blocks
163  int NumBlocks() const { return nBlocks; }
164 
165  //! Return a reference to block i,i.
167  { MFEM_VERIFY(ops[iblock], ""); return *ops[iblock]; }
168 
169  //! Return a reference to block i,i (const version).
170  const Operator & GetDiagonalBlock(int iblock) const
171  { MFEM_VERIFY(ops[iblock], ""); return *ops[iblock]; }
172 
173  //! Return the offsets for block starts
174  Array<int> & Offsets() { return offsets; }
175 
176  //! Read only access to the offsets for block starts
177  const Array<int> & Offsets() const { return offsets; }
178 
179  /// Operator application
180  virtual void Mult (const Vector & x, Vector & y) const;
181 
182  /// Action of the transpose operator
183  virtual void MultTranspose (const Vector & x, Vector & y) const;
184 
186 
187  //! Controls the ownership of the blocks: if nonzero,
188  //! BlockDiagonalPreconditioner will delete all blocks that are set
189  //! (non-NULL); the default value is zero.
191 
192 private:
193  //! Number of Blocks
194  int nBlocks;
195  //! Offsets for the starting position of each block
196  Array<int> offsets;
197  //! 1D array that stores each block of the operator.
198  Array<Operator *> ops;
199  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
200  //! methods.
201  mutable BlockVector xblock;
202  mutable BlockVector yblock;
203 };
204 
205 //! @class BlockLowerTriangularPreconditioner
206 /**
207  * \brief A class to handle Block lower triangular preconditioners in a
208  * matrix-free implementation.
209  *
210  * Usage:
211  * - Use the constructors to define the block structure
212  * - Use SetBlock() to fill the BlockOperator
213  * - Diagonal blocks of the preconditioner should approximate the inverses of
214  * the diagonal block of the matrix
215  * - Off-diagonal blocks of the preconditioner should match/approximate those of
216  * the original matrix
217  * - Use the method Mult() and MultTranspose() to apply the operator to a vector.
218  *
219  * If a diagonal block is not set, it is assumed to be an identity block, if an
220  * off-diagonal block is not set, it is assumed to be a zero block.
221  *
222  */
224 {
225 public:
226  //! Constructor for BlockLowerTriangularPreconditioners with the same
227  //! block-structure for rows and columns.
228  /**
229  * @param offsets Offsets that mark the start of each row/column block
230  * (size nBlocks+1).
231  *
232  * @note BlockLowerTriangularPreconditioner will not own/copy the data
233  * contained in @a offsets.
234  */
236 
237  //! Add block op in the block-entry (iblock, iblock).
238  /**
239  * @param iblock The block will be inserted in location (iblock, iblock).
240  * @param op The Operator to be inserted.
241  */
242  void SetDiagonalBlock(int iblock, Operator *op);
243  //! Add a block opt in the block-entry (iblock, jblock).
244  /**
245  * @param iRow, iCol The block will be inserted in location (iRow, iCol).
246  * @param op The Operator to be inserted.
247  */
248  void SetBlock(int iRow, int iCol, Operator *op);
249  //! This method is present since required by the abstract base class Solver
250  virtual void SetOperator(const Operator &op) { }
251 
252  //! Return the number of blocks
253  int NumBlocks() const { return nBlocks; }
254 
255  //! Return a reference to block i,j.
256  Operator & GetBlock(int iblock, int jblock)
257  { MFEM_VERIFY(ops(iblock,jblock), ""); return *ops(iblock,jblock); }
258 
259  //! Return the offsets for block starts
260  Array<int> & Offsets() { return offsets; }
261 
262  /// Operator application
263  virtual void Mult (const Vector & x, Vector & y) const;
264 
265  /// Action of the transpose operator
266  virtual void MultTranspose (const Vector & x, Vector & y) const;
267 
269 
270  //! Controls the ownership of the blocks: if nonzero,
271  //! BlockLowerTriangularPreconditioner will delete all blocks that are set
272  //! (non-NULL); the default value is zero.
274 
275 private:
276  //! Number of block rows/columns
277  int nBlocks;
278  //! Offsets for the starting position of each block
279  Array<int> offsets;
280  //! 2D array that stores each block of the operator.
282 
283  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
284  //! methods.
285  mutable BlockVector xblock;
286  mutable BlockVector yblock;
287  mutable Vector tmp;
288  mutable Vector tmp2;
289 };
290 
291 }
292 
293 #endif /* MFEM_BLOCKOPERATOR */
Array< int > & Offsets()
Return the offsets for block starts.
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
Array< int > & RowOffsets()
Return the row offsets for block starts.
Array< int > & Offsets()
Return the offsets for block starts.
BlockOperator & operator=(const BlockOperator &)=delete
Copy assignment is not supported.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
virtual void SetOperator(const Operator &op)
This method is present since required by the abstract base class Solver.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
int NumBlocks() const
Return the number of blocks.
A class to handle Block lower triangular preconditioners in a matrix-free implementation.
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
void SetBlock(int iRow, int iCol, Operator *op)
Add a block opt in the block-entry (iblock, jblock).
Operator & GetBlock(int iblock, int jblock)
Return a reference to block i,j.
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
BlockDiagonalPreconditioner(const Array< int > &offsets)
Constructor that specifies the block structure.
Dynamic 2D array using row-major layout.
Definition: array.hpp:351
BlockOperator(const Array< int > &offsets)
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
void SetDiagonalBlock(int iblock, Operator *op, double c=1.0)
Add block op in the block-entry (iblock, iblock).
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void SetBlockCoef(int i, int j, double c)
Set the coefficient for block i,j.
virtual void SetOperator(const Operator &op)
This method is present since required by the abstract base class Solver.
const Array< int > & RowOffsets() const
Read only access to the row offsets for block starts.
const Array< int > & ColOffsets() const
Read only access to the columns offsets for block starts.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void SetDiagonalBlock(int iblock, Operator *op)
Add block op in the block-entry (iblock, iblock).
int NumBlocks() const
Return the number of blocks.
const Operator & GetBlock(int i, int j) const
Return a reference to block i,j (const version)
const Array< int > & Offsets() const
Read only access to the offsets for block starts.
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
double GetBlockCoef(int i, int j) const
Return the coefficient for block i,j.
int NumRowBlocks() const
Return the number of row blocks.
Vector data type.
Definition: vector.hpp:60
BlockLowerTriangularPreconditioner(const Array< int > &offsets)
Base class for solvers.
Definition: operator.hpp:651
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
const Operator & GetDiagonalBlock(int iblock) const
Return a reference to block i,i (const version).
int NumColBlocks() const
Return the number of column blocks.
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Array< int > & ColOffsets()
Return the columns offsets for block starts.