MFEM  v4.5.2
Finite element discretization library
blockoperator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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  virtual Type GetType() const { return MFEM_Block_Operator; }
117 
118 private:
119  //! Number of block rows
120  int nRowBlocks;
121  //! Number of block columns
122  int nColBlocks;
123  //! Row offsets for the starting position of each block
124  Array<int> row_offsets;
125  //! Column offsets for the starting position of each block
126  Array<int> col_offsets;
127  //! 2D array that stores each block of the operator.
129  //! 2D array that stores a coefficient for each block of the operator.
130  Array2D<double> coef;
131 
132  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose methods.
133  mutable BlockVector xblock;
134  mutable BlockVector yblock;
135  mutable Vector tmp;
136 };
137 
138 //! @class BlockDiagonalPreconditioner
139 /**
140  * \brief A class to handle Block diagonal preconditioners in a matrix-free implementation.
141  *
142  * Usage:
143  * - Use the constructors to define the block structure
144  * - Use SetDiagonalBlock to fill the BlockOperator
145  * - Use the method Mult and MultTranspose to apply the operator to a vector.
146  *
147  * If a block is not set, it is assumed to be an identity block.
148  *
149  */
151 {
152 public:
153  //! Constructor that specifies the block structure
154  BlockDiagonalPreconditioner(const Array<int> & offsets);
155  //! Add a square block op in the block-entry (iblock, iblock).
156  /**
157  * iblock: The block will be inserted in location (iblock, iblock).
158  * op: the Operator to be inserted.
159  */
160  void SetDiagonalBlock(int iblock, Operator *op);
161  //! This method is present since required by the abstract base class Solver
162  virtual void SetOperator(const Operator &op) { }
163 
164  //! Return the number of blocks
165  int NumBlocks() const { return nBlocks; }
166 
167  //! Return a reference to block i,i.
169  { MFEM_VERIFY(ops[iblock], ""); return *ops[iblock]; }
170 
171  //! Return a reference to block i,i (const version).
172  const Operator & GetDiagonalBlock(int iblock) const
173  { MFEM_VERIFY(ops[iblock], ""); return *ops[iblock]; }
174 
175  //! Return the offsets for block starts
176  Array<int> & Offsets() { return offsets; }
177 
178  //! Read only access to the offsets for block starts
179  const Array<int> & Offsets() const { return offsets; }
180 
181  /// Operator application
182  virtual void Mult (const Vector & x, Vector & y) const;
183 
184  /// Action of the transpose operator
185  virtual void MultTranspose (const Vector & x, Vector & y) const;
186 
188 
189  //! Controls the ownership of the blocks: if nonzero,
190  //! BlockDiagonalPreconditioner will delete all blocks that are set
191  //! (non-NULL); the default value is zero.
193 
194 private:
195  //! Number of Blocks
196  int nBlocks;
197  //! Offsets for the starting position of each block
198  Array<int> offsets;
199  //! 1D array that stores each block of the operator.
200  Array<Operator *> ops;
201  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
202  //! methods.
203  mutable BlockVector xblock;
204  mutable BlockVector yblock;
205 };
206 
207 //! @class BlockLowerTriangularPreconditioner
208 /**
209  * \brief A class to handle Block lower triangular preconditioners in a
210  * matrix-free implementation.
211  *
212  * Usage:
213  * - Use the constructors to define the block structure
214  * - Use SetBlock() to fill the BlockOperator
215  * - Diagonal blocks of the preconditioner should approximate the inverses of
216  * the diagonal block of the matrix
217  * - Off-diagonal blocks of the preconditioner should match/approximate those of
218  * the original matrix
219  * - Use the method Mult() and MultTranspose() to apply the operator to a vector.
220  *
221  * If a diagonal block is not set, it is assumed to be an identity block, if an
222  * off-diagonal block is not set, it is assumed to be a zero block.
223  *
224  */
226 {
227 public:
228  //! Constructor for BlockLowerTriangularPreconditioners with the same
229  //! block-structure for rows and columns.
230  /**
231  * @param offsets Offsets that mark the start of each row/column block
232  * (size nBlocks+1).
233  *
234  * @note BlockLowerTriangularPreconditioner will not own/copy the data
235  * contained in @a offsets.
236  */
238 
239  //! Add block op in the block-entry (iblock, iblock).
240  /**
241  * @param iblock The block will be inserted in location (iblock, iblock).
242  * @param op The Operator to be inserted.
243  */
244  void SetDiagonalBlock(int iblock, Operator *op);
245  //! Add a block opt in the block-entry (iblock, jblock).
246  /**
247  * @param iRow, iCol The block will be inserted in location (iRow, iCol).
248  * @param op The Operator to be inserted.
249  */
250  void SetBlock(int iRow, int iCol, Operator *op);
251  //! This method is present since required by the abstract base class Solver
252  virtual void SetOperator(const Operator &op) { }
253 
254  //! Return the number of blocks
255  int NumBlocks() const { return nBlocks; }
256 
257  //! Return a reference to block i,j.
258  Operator & GetBlock(int iblock, int jblock)
259  { MFEM_VERIFY(ops(iblock,jblock), ""); return *ops(iblock,jblock); }
260 
261  //! Return the offsets for block starts
262  Array<int> & Offsets() { return offsets; }
263 
264  /// Operator application
265  virtual void Mult (const Vector & x, Vector & y) const;
266 
267  /// Action of the transpose operator
268  virtual void MultTranspose (const Vector & x, Vector & y) const;
269 
271 
272  //! Controls the ownership of the blocks: if nonzero,
273  //! BlockLowerTriangularPreconditioner will delete all blocks that are set
274  //! (non-NULL); the default value is zero.
276 
277 private:
278  //! Number of block rows/columns
279  int nBlocks;
280  //! Offsets for the starting position of each block
281  Array<int> offsets;
282  //! 2D array that stores each block of the operator.
284 
285  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
286  //! methods.
287  mutable BlockVector xblock;
288  mutable BlockVector yblock;
289  mutable Vector tmp;
290  mutable Vector tmp2;
291 };
292 
293 }
294 
295 #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.
double GetBlockCoef(int i, int j) const
Return the coefficient for block i,j.
int NumBlocks() const
Return the number of blocks.
Array< int > & Offsets()
Return the offsets for block starts.
ID for the base class BlockOperator.
Definition: operator.hpp:299
int NumRowBlocks() const
Return the number of row blocks.
BlockOperator & operator=(const BlockOperator &)=delete
Copy assignment is not supported.
const Operator & GetBlock(int i, int j) const
Return a reference to block i,j (const version)
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual void SetOperator(const Operator &op)
This method is present since required by the abstract base class Solver.
int NumColBlocks() const
Return the number of column blocks.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
A class to handle Block lower triangular preconditioners in a matrix-free implementation.
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
int NumBlocks() const
Return the number of blocks.
const Array< int > & RowOffsets() const
Read only access to the row offsets for block starts.
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.
BlockDiagonalPreconditioner(const Array< int > &offsets)
Constructor that specifies the block structure.
const Operator & GetDiagonalBlock(int iblock) const
Return a reference to block i,i (const version).
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
BlockOperator(const Array< int > &offsets)
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator.
void SetDiagonalBlock(int iblock, Operator *op, double c=1.0)
Add block op in the block-entry (iblock, iblock).
void SetBlockCoef(int i, int j, double c)
Set the coefficient for block i,j.
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
virtual Type GetType() const
virtual void SetOperator(const Operator &op)
This method is present since required by the abstract base class Solver.
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
void SetDiagonalBlock(int iblock, Operator *op)
Add block op in the block-entry (iblock, iblock).
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.
Vector data type.
Definition: vector.hpp:60
BlockLowerTriangularPreconditioner(const Array< int > &offsets)
const Array< int > & Offsets() const
Read only access to the offsets for block starts.
Base class for solvers.
Definition: operator.hpp:682
Abstract operator.
Definition: operator.hpp:24
A class to handle Block systems in a matrix-free implementation.
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).
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Array< int > & ColOffsets()
Return the columns offsets for block starts.