MFEM  v3.4
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
blockmatrix.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_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  if (row_offsets[iblock+1] > iglobal)
175  {
176  break;
177  }
178 
179  iloc = iglobal - row_offsets[iblock];
180 }
181 
182 inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock,
183  int & jloc) const
184 {
185  if (jglobal > col_offsets[nColBlocks])
186  {
187  mfem_error("BlockMatrix::findGlobalCol");
188  }
189 
190  for (jblock = 0; jblock < nColBlocks; ++jblock)
191  if (col_offsets[jblock+1] > jglobal)
192  {
193  break;
194  }
195 
196  jloc = jglobal - col_offsets[jblock];
197 }
198 
199 }
200 
201 #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:478
Abstract data type for sparse matrices.
Definition: matrix.hpp:77
Abstract data type for matrix inverse.
Definition: matrix.hpp:66
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:60
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:38
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:146
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:418
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:289
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:81
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:64
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...