MFEM  v4.6.0
Finite element discretization library
blockmatrix.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_BLOCKMATRIX
13 #define MFEM_BLOCKMATRIX
14 
15 #include "../config/config.hpp"
16 #include "../general/array.hpp"
17 #include "../general/globals.hpp"
18 #include "vector.hpp"
19 #include "sparsemat.hpp"
20 
21 namespace mfem
22 {
23 
25 {
26 public:
27  //! Constructor for square block matrices
28  /**
29  @param offsets offsets that mark the start of each row/column block (size
30  nRowBlocks+1).
31  @note BlockMatrix will not own/copy the data contained in offsets.
32  */
33  BlockMatrix(const Array<int> & offsets);
34  //! Constructor for rectangular block matrices
35  /**
36  @param row_offsets offsets that mark the start of each row block (size
37  nRowBlocks+1).
38  @param col_offsets offsets that mark the start of each column block (size
39  nColBlocks+1).
40  @note BlockMatrix will not own/copy the data contained in offsets.
41  */
42  BlockMatrix(const Array<int> & row_offsets, const Array<int> & col_offsets);
43  //! Set A(i,j) = mat
44  void SetBlock(int i, int j, SparseMatrix * mat);
45  //! Return the number of row blocks
46  int NumRowBlocks() const {return nRowBlocks; }
47  //! Return the number of column blocks
48  int NumColBlocks() const {return nColBlocks; }
49  //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL
50  SparseMatrix & GetBlock(int i, int j);
51  //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j)
52  //! == NULL. (const version)
53  const SparseMatrix & GetBlock(int i, int j) const;
54  //! Check if block (i,j) is a zero block.
55  int IsZeroBlock(int i, int j) const {return (Aij(i,j)==NULL) ? 1 : 0; }
56  //! Return the row offsets for block starts
57  Array<int> & RowOffsets() { return row_offsets; }
58  //! Return the columns offsets for block starts
59  Array<int> & ColOffsets() { return col_offsets; }
60  //! Return the row offsets for block starts (const version)
61  const Array<int> & RowOffsets() const { return row_offsets; }
62  //! Return the row offsets for block starts (const version)
63  const Array<int> & ColOffsets() const { return col_offsets; }
64  //! Return the number of non zeros in row i
65  int RowSize(const int i) const;
66 
67  /// Eliminate the row and column @a rc from the matrix.
68  /** Eliminates the column and row @a rc, replacing the element (rc,rc) with
69  1.0. Assumes that element (i,rc) is assembled if and only if the element
70  (rc,i) is assembled. If @a dpolicy is specified, the element (rc,rc) is
71  treated according to that policy. */
72  void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE);
73 
74  /** @brief Eliminate the rows and columns corresponding to the entries
75  in @a vdofs + save the eliminated entries into
76  @a Ae so that (*this) + Ae is equal to the original matrix. */
77  void EliminateRowCols(const Array<int> & vdofs, BlockMatrix *Ae,
78  DiagonalPolicy dpolicy = DIAG_ONE);
79 
80  //! Symmetric elimination of the marked degree of freedom.
81  /**
82  @param ess_bc_dofs marker of the degree of freedom to be eliminated
83  dof i is eliminated if @a ess_bc_dofs[i] = 1.
84  @param sol vector that stores the values of the degree of freedom
85  that need to be eliminated
86  @param rhs vector that stores the rhs of the system.
87  */
88  void EliminateRowCol(Array<int> & ess_bc_dofs, Vector & sol, Vector & rhs);
89 
90  /// Finalize all the submatrices
91  virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
92  /// A slightly more general version of the Finalize(int) method.
93  void Finalize(int skip_zeros, bool fix_empty_rows);
94 
95  //! Returns a monolithic CSR matrix that represents this operator.
97  //! Export the monolithic matrix to file.
98  virtual void PrintMatlab(std::ostream & os = mfem::out) const;
99 
100  /// @name Matrix interface
101  ///@{
102 
103  /// Returns reference to a_{ij}.
104  virtual double& Elem (int i, int j);
105  /// Returns constant reference to a_{ij}.
106  virtual const double& Elem (int i, int j) const;
107  /// Returns a pointer to (approximation) of the matrix inverse.
108  virtual MatrixInverse * Inverse() const
109  {
110  mfem_error("BlockMatrix::Inverse not implemented \n");
111  return static_cast<MatrixInverse*>(NULL);
112  }
113  ///@}
114 
115  ///@name AbstractSparseMatrix interface
116  ///@{
117 
118  //! Returns the total number of non zeros in the matrix.
119  virtual int NumNonZeroElems() const;
120  /// Gets the columns indexes and values for row *row*.
121  /** The return value is always 0 since @a cols and @a srow are copies of the
122  values in the matrix. */
123  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
124  /** @brief If the matrix is square, this method will place 1 on the diagonal
125  (i,i) if row i has "almost" zero l1-norm.
126 
127  If entry (i,i) does not belong to the sparsity pattern of A, then a error
128  will occur. */
129  virtual void EliminateZeroRows(const double threshold = 1e-12);
130 
131  /// Matrix-Vector Multiplication y = A*x
132  virtual void Mult(const Vector & x, Vector & y) const;
133  /// Matrix-Vector Multiplication y = y + val*A*x
134  virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const;
135  /// MatrixTranspose-Vector Multiplication y = A'*x
136  virtual void MultTranspose(const Vector & x, Vector & y) const;
137  /// MatrixTranspose-Vector Multiplication y = y + val*A'*x
138  virtual void AddMultTranspose(const Vector & x, Vector & y,
139  const double val = 1.) const;
140  ///@}
141 
142  /** @brief Partial matrix vector multiplication of (*this) with @a x
143  involving only the rows given by @a rows. The result is given in @a y */
144  void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
145  /** @brief Partial matrix vector multiplication of (*this) with @a x
146  involving only the rows given by @a rows. The result is multiplied by
147  @a a and added to @a y */
148  void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
149  const double a=1.0) const;
150 
151  //! Destructor
152  virtual ~BlockMatrix();
153  //! If owns_blocks the SparseMatrix objects Aij will be deallocated.
155 
156  virtual Type GetType() const { return MFEM_Block_Matrix; }
157 
158 private:
159  //! Given a global row iglobal finds to which row iloc in block iblock belongs to.
160  inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const;
161  //! Given a global column jglobal finds to which column jloc in block jblock belongs to.
162  inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const;
163 
164  //! Number of row blocks
165  int nRowBlocks;
166  //! Number of columns blocks
167  int nColBlocks;
168  //! row offsets for each block start (length nRowBlocks+1).
169  Array<int> row_offsets;
170  //! column offsets for each block start (length nColBlocks+1).
171  Array<int> col_offsets;
172  //! 2D array that stores each block of the BlockMatrix. Aij(iblock, jblock)
173  //! == NULL if block (iblock, jblock) is all zeros.
175 };
176 
177 //! Transpose a BlockMatrix: result = A'
178 BlockMatrix * Transpose(const BlockMatrix & A);
179 //! Multiply BlockMatrix matrices: result = A*B
180 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B);
181 
182 inline void BlockMatrix::findGlobalRow(int iglobal, int & iblock,
183  int & iloc) const
184 {
185  if (iglobal > row_offsets[nRowBlocks])
186  {
187  mfem_error("BlockMatrix::findGlobalRow");
188  }
189 
190  for (iblock = 0; iblock < nRowBlocks; ++iblock)
191  {
192  if (row_offsets[iblock+1] > iglobal) { break; }
193  }
194 
195  iloc = iglobal - row_offsets[iblock];
196 }
197 
198 inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock,
199  int & jloc) const
200 {
201  if (jglobal > col_offsets[nColBlocks])
202  {
203  mfem_error("BlockMatrix::findGlobalCol");
204  }
205 
206  for (jblock = 0; jblock < nColBlocks; ++jblock)
207  {
208  if (col_offsets[jblock+1] > jglobal) { break; }
209  }
210 
211  jloc = jglobal - col_offsets[jblock];
212 }
213 
214 }
215 
216 #endif /* MFEM_BLOCKMATRIX */
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows...
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse.
const Array< int > & RowOffsets() const
Return the row offsets for block starts (const version)
Definition: blockmatrix.hpp:61
int owns_blocks
If owns_blocks the SparseMatrix objects Aij will be deallocated.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const double a=1.0) const
Partial matrix vector multiplication of (*this) with x involving only the rows given by rows...
Abstract data type for sparse matrices.
Definition: matrix.hpp:73
virtual void PrintMatlab(std::ostream &os=mfem::out) const
Export the monolithic matrix to file.
Abstract data type for matrix inverse.
Definition: matrix.hpp:62
ID for class BlockMatrix.
Definition: operator.hpp:298
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:23
const Array< int > & ColOffsets() const
Return the row offsets for block starts (const version)
Definition: blockmatrix.hpp:63
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:61
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:46
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual void Finalize(int skip_zeros=1)
Finalize all the submatrices.
Definition: blockmatrix.hpp:91
Data type sparse matrix.
Definition: sparsemat.hpp:50
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A&#39;*x.
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:59
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:48
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:57
Set the diagonal value to one.
Definition: operator.hpp:50
virtual Type GetType() const
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:49
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Gets the columns indexes and values for row row.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
int RowSize(const int i) const
Return the number of non zeros in row i.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
void EliminateRowCol(int rc, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the row and column rc from the matrix.
double a
Definition: lissajous.cpp:41
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Definition: blockmatrix.cpp:82
virtual void AddMult(const Vector &x, Vector &y, const double val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
Vector data type.
Definition: vector.hpp:58
SparseMatrix * CreateMonolithic() const
Returns a monolithic CSR matrix that represents this operator.
virtual void AddMultTranspose(const Vector &x, Vector &y, const double val=1.) const
MatrixTranspose-Vector Multiplication y = y + val*A&#39;*x.
virtual void EliminateZeroRows(const double threshold=1e-12)
If the matrix is square, this method will place 1 on the diagonal (i,i) if row i has "almost" zero l1...
void EliminateRowCols(const Array< int > &vdofs, BlockMatrix *Ae, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the rows and columns corresponding to the entries in vdofs + save the eliminated entries in...
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockmatrix.hpp:55