MFEM  v4.1.0 Finite element discretization library
blockmatrix.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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  //! Symmetric elimination of the marked degree of freedom.
75  /**
76  @param ess_bc_dofs marker of the degree of freedom to be eliminated
77  dof i is eliminated if @a ess_bc_dofs[i] = 1.
78  @param sol vector that stores the values of the degree of freedom
79  that need to be eliminated
80  @param rhs vector that stores the rhs of the system.
81  */
82  void EliminateRowCol(Array<int> & ess_bc_dofs, Vector & sol, Vector & rhs);
83
84  /// Finalize all the submatrices
85  virtual void Finalize(int skip_zeros = 1) { Finalize(skip_zeros, false); }
86  /// A slightly more general version of the Finalize(int) method.
87  void Finalize(int skip_zeros, bool fix_empty_rows);
88
89  //! Returns a monolithic CSR matrix that represents this operator.
91  //! Export the monolithic matrix to file.
92  void PrintMatlab(std::ostream & os = mfem::out) const;
93
94  /// @name Matrix interface
95  ///@{
96
97  /// Returns reference to a_{ij}.
98  virtual double& Elem (int i, int j);
99  /// Returns constant reference to a_{ij}.
100  virtual const double& Elem (int i, int j) const;
101  /// Returns a pointer to (approximation) of the matrix inverse.
102  virtual MatrixInverse * Inverse() const
103  {
104  mfem_error("BlockMatrix::Inverse not implemented \n");
105  return static_cast<MatrixInverse*>(NULL);
106  }
107  ///@}
108
109  ///@name AbstractSparseMatrix interface
110  ///@{
111
112  //! Returns the total number of non zeros in the matrix.
113  virtual int NumNonZeroElems() const;
114  /// Gets the columns indexes and values for row *row*.
115  /** The return value is always 0 since @a cols and @a srow are copies of the
116  values in the matrix. */
117  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
118  /** @brief If the matrix is square, this method will place 1 on the diagonal
119  (i,i) if row i has "almost" zero l1-norm.
120
121  If entry (i,i) does not belong to the sparsity pattern of A, then a error
122  will occur. */
123  virtual void EliminateZeroRows(const double threshold = 1e-12);
124
125  /// Matrix-Vector Multiplication y = A*x
126  virtual void Mult(const Vector & x, Vector & y) const;
127  /// Matrix-Vector Multiplication y = y + val*A*x
128  virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const;
129  /// MatrixTranspose-Vector Multiplication y = A'*x
130  virtual void MultTranspose(const Vector & x, Vector & y) const;
131  /// MatrixTranspose-Vector Multiplication y = y + val*A'*x
132  virtual void AddMultTranspose(const Vector & x, Vector & y,
133  const double val = 1.) const;
134  ///@}
135
136  //! Destructor
137  virtual ~BlockMatrix();
138  //! If owns_blocks the SparseMatrix objects Aij will be deallocated.
140
141 private:
142  //! Given a global row iglobal finds to which row iloc in block iblock belongs to.
143  inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const;
144  //! Given a global column jglobal finds to which column jloc in block jblock belongs to.
145  inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const;
146
147  //! Number of row blocks
148  int nRowBlocks;
149  //! Number of columns blocks
150  int nColBlocks;
151  //! row offsets for each block start (length nRowBlocks+1).
152  Array<int> row_offsets;
153  //! column offsets for each block start (length nColBlocks+1).
154  Array<int> col_offsets;
155  //! 2D array that stores each block of the BlockMatrix. Aij(iblock, jblock)
156  //! == NULL if block (iblock, jblock) is all zeros.
158 };
159
160 //! Transpose a BlockMatrix: result = A'
161 BlockMatrix * Transpose(const BlockMatrix & A);
162 //! Multiply BlockMatrix matrices: result = A*B
163 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B);
164
165 inline void BlockMatrix::findGlobalRow(int iglobal, int & iblock,
166  int & iloc) const
167 {
168  if (iglobal > row_offsets[nRowBlocks])
169  {
170  mfem_error("BlockMatrix::findGlobalRow");
171  }
172
173  for (iblock = 0; iblock < nRowBlocks; ++iblock)
174  {
175  if (row_offsets[iblock+1] > iglobal) { break; }
176  }
177
178  iloc = iglobal - row_offsets[iblock];
179 }
180
181 inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock,
182  int & jloc) const
183 {
184  if (jglobal > col_offsets[nColBlocks])
185  {
186  mfem_error("BlockMatrix::findGlobalCol");
187  }
188
189  for (jblock = 0; jblock < nColBlocks; ++jblock)
190  {
191  if (col_offsets[jblock+1] > jglobal) { break; }
192  }
193
194  jloc = jglobal - col_offsets[jblock];
195 }
196
197 }
198
199 #endif /* MFEM_BLOCKMATRIX */
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse.
SparseMatrix * CreateMonolithic() const
Returns a monolithic CSR matrix that represents this operator.
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:476
Abstract data type for sparse matrices.
Definition: matrix.hpp:80
Abstract data type for matrix inverse.
Definition: matrix.hpp:69
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:22
virtual void AddMultTranspose(const Vector &x, Vector &y, const double val=1.) const
MatrixTranspose-Vector Multiplication y = y + val*A&#39;*x.
virtual int GetRow(const int row, Array< int > &cols, Vector &srow) const
Gets the columns indexes and values for row row.
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:62
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:85
Data type sparse matrix.
Definition: sparsemat.hpp:40
int IsZeroBlock(int i, int j) const
Check if block (i,j) is a zero block.
Definition: blockmatrix.hpp:55
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:153
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:46
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:59
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:414
void PrintMatlab(std::ostream &os=mfem::out) const
Export the monolithic matrix to file.
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:57
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:50
Set the diagonal value to one.
Definition: matrix.hpp:35
Dynamic 2D array using row-major layout.
Definition: array.hpp:316
virtual void MultTranspose(const Vector &x, Vector &y) const
MatrixTranspose-Vector Multiplication y = A&#39;*x.
virtual void Mult(const Vector &x, Vector &y) const
Matrix-Vector Multiplication y = A*x.
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:48
void EliminateRowCol(int rc, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the row and column rc from the matrix.
virtual void AddMult(const Vector &x, Vector &y, const double val=1.) const
Matrix-Vector Multiplication y = y + val*A*x.
const Array< int > & ColOffsets() const
Return the row offsets for block starts (const version)
Definition: blockmatrix.hpp:63
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:83
Vector data type.
Definition: vector.hpp:48
const Array< int > & RowOffsets() const
Return the row offsets for block starts (const version)
Definition: blockmatrix.hpp:61
int RowSize(const int i) const
Return the number of non zeros in row i.
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
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 &quot;almost&quot; zero l1...