MFEM  v3.3 Finite element discretization library
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
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 "vector.hpp"
18 #include "sparsemat.hpp"
19
20 namespace mfem
21 {
22
24 {
25 public:
26  //! Constructor for square block matrices
27  /**
28  * offsets: offsets that mark the start of each row/column block (size nRowBlocks+1).
29  * Note: BlockMatrix will not own/copy the data contained in offsets.
30  */
31  BlockMatrix(const Array<int> & offsets);
32  //! Constructor for rectangular block matrices
33  /**
34  * row_offsets: offsets that mark the start of each row block (size nRowBlocks+1).
35  * col_offsets: offsets that mark the start of each column block (size nColBlocks+1).
36  * Note: BlockMatrix will not own/copy the data contained in offsets.
37  */
38  BlockMatrix(const Array<int> & row_offsets, const Array<int> & col_offsets);
39  //! Set A(i,j) = mat
40  void SetBlock(int i, int j, SparseMatrix * mat);
41  //! Return the number of row blocks
42  int NumRowBlocks() const {return nRowBlocks; }
43  //! Return the number of column blocks
44  int NumColBlocks() const {return nColBlocks; }
45  //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL
46  SparseMatrix & GetBlock(int i, int j);
47  //! Return a reference to block (i,j). Reference may be invalid if Aij(i,j)
48  //! == NULL. (const version)
49  const SparseMatrix & GetBlock(int i, int j) const;
50  //! Check if block (i,j) is a zero block.
51  int IsZeroBlock(int i, int j) const {return (Aij(i,j)==NULL) ? 1 : 0; }
52  //! Return the row offsets for block starts
53  Array<int> & RowOffsets() { return row_offsets; }
54  //! Return the columns offsets for block starts
55  Array<int> & ColOffsets() { return col_offsets; }
56  //! Return the row offsets for block starts (const version)
57  const Array<int> & RowOffsets() const { return row_offsets; }
58  //! Return the row offsets for block starts (const version)
59  const Array<int> & ColOffsets() const { return col_offsets; }
60  //! Return the number of non zeros in row i
61  int RowSize(const int i) const;
62  //! Symmetric elimination of the marked degree of freedom.
63  /**
64  * ess_bc_dofs: marker of the degree of freedom to be eliminated
65  * dof i is eliminated if ess_bc_dofs[i] = 1.
66  * sol: vector that stores the values of the degree of freedom that need to
67  * be eliminated
68  * rhs: vector that stores the rhs of the system.
69  */
70  void EliminateRowCol(Array<int> & ess_bc_dofs, Vector & sol, Vector & rhs);
71  //! Returns a monolithic CSR matrix that represents this operator.
73  //! Export the monolithic matrix to file.
74  void PrintMatlab(std::ostream & os = std::cout) const;
75
76  //@name Matrix interface
77  //@{
78  /// Returns reference to a_{ij}.
79  virtual double& Elem (int i, int j);
80  /// Returns constant reference to a_{ij}.
81  virtual const double& Elem (int i, int j) const;
82  /// Returns a pointer to (approximation) of the matrix inverse.
83  virtual MatrixInverse * Inverse() const
84  {
85  mfem_error("BlockMatrix::Inverse not implemented \n");
86  return static_cast<MatrixInverse*>(NULL);
87  }
88  //@}
89
90  //@name AbstractSparseMatrix interface
91  //@{
92  //! Returns the total number of non zeros in the matrix.
93  virtual int NumNonZeroElems() const;
94  /// Gets the columns indexes and values for row *row*.
95  /// The return value is always 0 since cols and srow are copies of the values in the matrix.
96  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
97  //! If the matrix is square, it will place 1 on the diagonal (i,i) if row i
98  //! has "almost" zero l1-norm.
99  /**
100  * If entry (i,i) does not belong to the sparsity pattern of A, then a error will occur.
101  */
102  virtual void EliminateZeroRows();
103
104  /// Matrix-Vector Multiplication y = A*x
105  virtual void Mult(const Vector & x, Vector & y) const;
106  /// Matrix-Vector Multiplication y = y + val*A*x
107  virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const;
108  /// MatrixTranspose-Vector Multiplication y = A'*x
109  virtual void MultTranspose(const Vector & x, Vector & y) const;
110  /// MatrixTranspose-Vector Multiplication y = y + val*A'*x
111  virtual void AddMultTranspose(const Vector & x, Vector & y,
112  const double val = 1.) const;
113  //@}
114
115  //! Destructor
116  virtual ~BlockMatrix();
117  //! if owns_blocks the SparseMatrix objects Aij will be deallocated.
119
120 private:
121  //! Given a global row iglobal finds to which row iloc in block iblock belongs to.
122  inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const;
123  //! Given a global column jglobal finds to which column jloc in block jblock belongs to.
124  inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const;
125
126  //! Number of row blocks
127  int nRowBlocks;
128  //! Number of columns blocks
129  int nColBlocks;
130  //! row offsets for each block start (length nRowBlocks+1).
131  Array<int> row_offsets;
132  //! column offsets for each block start (length nColBlocks+1).
133  Array<int> col_offsets;
134  //! 2D array that stores each block of the BlockMatrix. Aij(iblock, jblock)
135  //! == NULL if block (iblock, jblock) is all zeros.
137 };
138
139 //! Transpose a BlockMatrix: result = A'
140 BlockMatrix * Transpose(const BlockMatrix & A);
141 //! Multiply BlockMatrix matrices: result = A*B
142 BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B);
143
144 inline void BlockMatrix::findGlobalRow(int iglobal, int & iblock,
145  int & iloc) const
146 {
147  if (iglobal > row_offsets[nRowBlocks])
148  {
149  mfem_error("BlockMatrix::findGlobalRow");
150  }
151
152  for (iblock = 0; iblock < nRowBlocks; ++iblock)
153  if (row_offsets[iblock+1] > iglobal)
154  {
155  break;
156  }
157
158  iloc = iglobal - row_offsets[iblock];
159 }
160
161 inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock,
162  int & jloc) const
163 {
164  if (jglobal > col_offsets[nColBlocks])
165  {
166  mfem_error("BlockMatrix::findGlobalCol");
167  }
168
169  for (jblock = 0; jblock < nColBlocks; ++jblock)
170  if (col_offsets[jblock+1] > jglobal)
171  {
172  break;
173  }
174
175  jloc = jglobal - col_offsets[jblock];
176 }
177
178 }
179
180 #endif /* MFEM_BLOCKMATRIX */
void PrintMatlab(std::ostream &os=std::cout) const
Export the monolithic matrix to file.
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse.
Definition: blockmatrix.hpp:83
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:468
Abstract data type for sparse matrices.
Definition: matrix.hpp:69
Abstract data type for matrix inverse.
Definition: matrix.hpp:58
BlockMatrix(const Array< int > &offsets)
Constructor for square block matrices.
Definition: blockmatrix.cpp:21
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
void SetBlock(int i, int j, SparseMatrix *mat)
Set A(i,j) = mat.
Definition: blockmatrix.cpp:59
virtual double & Elem(int i, int j)
Returns reference to a_{ij}.
virtual void EliminateZeroRows()
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:51
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:42
Array< int > & ColOffsets()
Return the columns offsets for block starts.
Definition: blockmatrix.hpp:55
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
Definition: table.cpp:408
Array< int > & RowOffsets()
Return the row offsets for block starts.
Definition: blockmatrix.hpp:53
void EliminateRowCol(Array< int > &ess_bc_dofs, Vector &sol, Vector &rhs)
Symmetric elimination of the marked degree of freedom.
virtual int NumNonZeroElems() const
Returns the total number of non zeros in the matrix.
virtual ~BlockMatrix()
Destructor.
Definition: blockmatrix.cpp:49
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.
void mfem_error(const char *msg)
Definition: error.cpp:106
int NumColBlocks() const
Return the number of column blocks.
Definition: blockmatrix.hpp:44
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:59
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:80
Vector data type.
Definition: vector.hpp:36
const Array< int > & RowOffsets() const
Return the row offsets for block starts (const version)
Definition: blockmatrix.hpp:57
int RowSize(const int i) const
Return the number of non zeros in row i.