MFEM  v3.1
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 "vector.hpp"
18 #include "sparsemat.hpp"
19 
20 namespace mfem
21 {
22 
24 {
25 public:
27 
31  BlockMatrix(const Array<int> & offsets);
33 
38  BlockMatrix(const Array<int> & row_offsets, const Array<int> & col_offsets);
40  void SetBlock(int i, int j, SparseMatrix * mat);
42  int NumRowBlocks() const {return nRowBlocks; }
44  int NumColBlocks() const {return nColBlocks; }
46  SparseMatrix & GetBlock(int i, int j);
49  const SparseMatrix & GetBlock(int i, int j) const;
51  int IsZeroBlock(int i, int j) const {return (Aij(i,j)==NULL) ? 1 : 0; }
53  Array<int> & RowOffsets() { return row_offsets; }
55  Array<int> & ColOffsets() { return col_offsets; }
57  const Array<int> & RowOffsets() const { return row_offsets; }
59  const Array<int> & ColOffsets() const { return col_offsets; }
61  int RowSize(const int i) const;
63 
70  void EliminateRowCol(Array<int> & ess_bc_dofs, Vector & sol, Vector & rhs);
74  void PrintMatlab(std::ostream & os = std::cout) const;
75 
76  //@name Matrix interface
78  virtual double& Elem (int i, int j);
81  virtual const double& Elem (int i, int j) const;
83  virtual MatrixInverse * Inverse() const
84  {
85  mfem_error("BlockMatrix::Inverse not implemented \n");
86  return static_cast<MatrixInverse*>(NULL);
87  }
89 
90  //@name AbstractSparseMatrix interface
92  virtual int NumNonZeroElems() const;
96  virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;
99 
102  virtual void EliminateZeroRows();
103 
105  virtual void Mult(const Vector & x, Vector & y) const;
107  virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const;
109  virtual void MultTranspose(const Vector & x, Vector & y) const;
111  virtual void AddMultTranspose(const Vector & x, Vector & y,
112  const double val = 1.) const;
114 
116  virtual ~BlockMatrix();
119 
120 private:
122  inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const;
124  inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const;
125 
127  int nRowBlocks;
129  int nColBlocks;
131  Array<int> row_offsets;
133  Array<int> col_offsets;
137 };
138 
140 BlockMatrix * Transpose(const BlockMatrix & A);
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:445
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:385
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:23
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:33
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.