MFEM  v3.3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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  //! Add block op in the block-entry (iblock, iblock).
55  /**
56  * iblock: The block will be inserted in location (iblock, iblock).
57  * op: the Operator to be inserted.
58  * c: optional scalar multiple for this block.
59  */
60  void SetDiagonalBlock(int iblock, Operator *op, double c = 1.0);
61  //! Add a block op in the block-entry (iblock, jblock).
62  /**
63  * irow, icol: The block will be inserted in location (irow, icol).
64  * op: the Operator to be inserted.
65  * c: optional scalar multiple for this block.
66  */
67  void SetBlock(int iRow, int iCol, Operator *op, double c = 1.0);
68 
69  //! Return the number of row blocks
70  int NumRowBlocks() const { return nRowBlocks; }
71  //! Return the number of column blocks
72  int NumColBlocks() const { return nColBlocks; }
73 
74  //! Check if block (i,j) is a zero block
75  int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; }
76  //! Return a reference to block i,j
77  Operator & GetBlock(int i, int j)
78  { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
79  //! Return the coefficient for block i,j
80  double GetBlockCoef(int i, int j) const
81  { MFEM_VERIFY(op(i,j), ""); return coef(i,j); }
82  //! Set the coefficient for block i,j
83  void SetBlockCoef(int i, int j, double c)
84  { MFEM_VERIFY(op(i,j), ""); coef(i,j) = c; }
85 
86  //! Return the row offsets for block starts
87  Array<int> & RowOffsets() { return row_offsets; }
88  //! Return the columns offsets for block starts
89  Array<int> & ColOffsets() { return col_offsets; }
90 
91  /// Operator application
92  virtual void Mult (const Vector & x, Vector & y) const;
93 
94  /// Action of the transpose operator
95  virtual void MultTranspose (const Vector & x, Vector & y) const;
96 
98 
99  //! Controls the ownership of the blocks: if nonzero, BlockOperator will
100  //! delete all blocks that are set (non-NULL); the default value is zero.
102 
103 private:
104  //! Number of block rows
105  int nRowBlocks;
106  //! Number of block columns
107  int nColBlocks;
108  //! Row offsets for the starting position of each block
109  Array<int> row_offsets;
110  //! Column offsets for the starting position of each block
111  Array<int> col_offsets;
112  //! 2D array that stores each block of the operator.
114  //! 2D array that stores a coefficient for each block of the operator.
115  Array2D<double> coef;
116 
117  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose methods.
118  mutable BlockVector xblock;
119  mutable BlockVector yblock;
120  mutable Vector tmp;
121 };
122 
123 //! @class BlockDiagonalPreconditioner
124 /**
125  * \brief A class to handle Block diagonal preconditioners in a matrix-free implementation.
126  *
127  * Usage:
128  * - Use the constructors to define the block structure
129  * - Use SetDiagonalBlock to fill the BlockOperator
130  * - Use the method Mult and MultTranspose to apply the operator to a vector.
131  *
132  * If a block is not set, it is assumed to be an identity block.
133  *
134  */
136 {
137 public:
138  //! Constructor that specifies the block structure
139  BlockDiagonalPreconditioner(const Array<int> & offsets);
140  //! Add a square block op in the block-entry (iblock, iblock).
141  /**
142  * iblock: The block will be inserted in location (iblock, iblock).
143  * op: the Operator to be inserted.
144  */
145  void SetDiagonalBlock(int iblock, Operator *op);
146  //! This method is present since required by the abstract base class Solver
147  virtual void SetOperator(const Operator &op) { }
148 
149  //! Return the number of blocks
150  int NumBlocks() const { return nBlocks; }
151 
152  //! Return a reference to block i,i.
154  { MFEM_VERIFY(op[iblock], ""); return *op[iblock]; }
155 
156  //! Return the offsets for block starts
157  Array<int> & Offsets() { return offsets; }
158 
159  /// Operator application
160  virtual void Mult (const Vector & x, Vector & y) const;
161 
162  /// Action of the transpose operator
163  virtual void MultTranspose (const Vector & x, Vector & y) const;
164 
166 
167  //! Controls the ownership of the blocks: if nonzero,
168  //! BlockDiagonalPreconditioner will delete all blocks that are set
169  //! (non-NULL); the default value is zero.
171 
172 private:
173  //! Number of Blocks
174  int nBlocks;
175  //! Offsets for the starting position of each block
176  Array<int> offsets;
177  //! 1D array that stores each block of the operator.
179  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
180  //! methods.
181  mutable BlockVector xblock;
182  mutable BlockVector yblock;
183 };
184 
185 //! @class BlockLowerTriangularPreconditioner
186 /**
187  * \brief A class to handle Block lower triangular preconditioners in a
188  * matrix-free implementation.
189  *
190  * Usage:
191  * - Use the constructors to define the block structure
192  * - Use SetBlock() to fill the BlockOperator
193  * - Diagonal blocks of the preconditioner should approximate the inverses of
194  * the diagonal block of the matrix
195  * - Off-diagonal blocks of the preconditioner should match/approximate those of
196  * the original matrix
197  * - Use the method Mult() and MultTranspose() to apply the operator to a vector.
198  *
199  * If a diagonal block is not set, it is assumed to be an identity block, if an
200  * off-diagonal block is not set, it is assumed to be a zero block.
201  *
202  */
204 {
205 public:
206  //! Constructor for BlockLowerTriangularPreconditioners with the same
207  //! block-structure for rows and columns.
208  /**
209  * @param offsets Offsets that mark the start of each row/column block
210  * (size nBlocks+1).
211  *
212  * @note BlockLowerTriangularPreconditioner will not own/copy the data
213  * contained in @a offsets.
214  */
216 
217  //! Add block op in the block-entry (iblock, iblock).
218  /**
219  * @param iblock The block will be inserted in location (iblock, iblock).
220  * @param op The Operator to be inserted.
221  */
222  void SetDiagonalBlock(int iblock, Operator *op);
223  //! Add a block op in the block-entry (iblock, jblock).
224  /**
225  * @param iRow, iCol The block will be inserted in location (iRow, iCol).
226  * @param op The Operator to be inserted.
227  */
228  void SetBlock(int iRow, int iCol, Operator *op);
229  //! This method is present since required by the abstract base class Solver
230  virtual void SetOperator(const Operator &op) { }
231 
232  //! Return the number of blocks
233  int NumBlocks() const { return nBlocks; }
234 
235  //! Return a reference to block i,j.
236  Operator & GetBlock(int iblock, int jblock)
237  { MFEM_VERIFY(op(iblock,jblock), ""); return *op(iblock,jblock); }
238 
239  //! Return the offsets for block starts
240  Array<int> & Offsets() { return offsets; }
241 
242  /// Operator application
243  virtual void Mult (const Vector & x, Vector & y) const;
244 
245  /// Action of the transpose operator
246  virtual void MultTranspose (const Vector & x, Vector & y) const;
247 
249 
250  //! Controls the ownership of the blocks: if nonzero,
251  //! BlockLowerTriangularPreconditioner will delete all blocks that are set
252  //! (non-NULL); the default value is zero.
254 
255 private:
256  //! Number of block rows/columns
257  int nBlocks;
258  //! Offsets for the starting position of each block
259  Array<int> offsets;
260  //! 2D array that stores each block of the operator.
262 
263  //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
264  //! methods.
265  mutable BlockVector xblock;
266  mutable BlockVector yblock;
267  mutable Vector tmp;
268  mutable Vector tmp2;
269 };
270 
271 }
272 #endif /* MFEM_BLOCKOPERATOR */
Array< int > & Offsets()
Return the offsets for block starts.
Array< int > & RowOffsets()
Return the row offsets for block starts.
Array< int > & Offsets()
Return the offsets for block starts.
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 op 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.
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.
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.
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:41
BlockLowerTriangularPreconditioner(const Array< int > &offsets)
Base class for solvers.
Definition: operator.hpp:259
Abstract operator.
Definition: operator.hpp:21
A class to handle Block systems in a matrix-free implementation.
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.