blockmatrix.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.googlecode.com.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_BLOCKMATRIX
#define MFEM_BLOCKMATRIX

#include "../config/config.hpp"
#include "../general/array.hpp"
#include "vector.hpp"
#include "sparsemat.hpp"

namespace mfem
{

class BlockMatrix : public AbstractSparseMatrix
{
public:

   BlockMatrix(const Array<int> & offsets);

   BlockMatrix(const Array<int> & row_offsets, const Array<int> & col_offsets);
   void SetBlock(int i, int j, SparseMatrix * mat);
   int NumRowBlocks() const {return nRowBlocks; }
   int NumColBlocks() const {return nColBlocks; }
   SparseMatrix & GetBlock(int i, int j);
   const SparseMatrix & GetBlock(int i, int j) const;
   int IsZeroBlock(int i, int j) const {return (Aij(i,j)==NULL) ? 1 : 0; }
   Array<int> & RowOffsets() { return row_offsets; }
   Array<int> & ColOffsets() { return col_offsets; }
   const Array<int> & RowOffsets() const { return row_offsets; }
   const Array<int> & ColOffsets() const { return col_offsets; }
   int RowSize(const int i) const;

   void EliminateRowCol(Array<int> & ess_bc_dofs, Vector & sol, Vector & rhs);
   SparseMatrix * CreateMonolithic() const;
   void PrintMatlab(std::ostream & os = std::cout) const;

   //@name Matrix interface
   virtual double& Elem (int i, int j);
   virtual const double& Elem (int i, int j) const;
   virtual MatrixInverse * Inverse() const
   {
      mfem_error("BlockMatrix::Inverse not implemented \n");
      return static_cast<MatrixInverse*>(NULL);
   }

   //@name AbstractSparseMatrix interface
   virtual int NumNonZeroElems() const;
   virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const;

   virtual void EliminateZeroRows();

   virtual void Mult(const Vector & x, Vector & y) const;
   virtual void AddMult(const Vector & x, Vector & y, const double val = 1.) const;
   virtual void MultTranspose(const Vector & x, Vector & y) const;
   virtual void AddMultTranspose(const Vector & x, Vector & y, const double val = 1.) const;

   virtual ~BlockMatrix();
   int owns_blocks;

private:
   inline void findGlobalRow(int iglobal, int & iblock, int & iloc) const;
   inline void findGlobalCol(int jglobal, int & jblock, int & jloc) const;

   int nRowBlocks;
   int nColBlocks;
   Array<int> row_offsets;
   Array<int> col_offsets;
   Array2D<SparseMatrix *> Aij;
};

BlockMatrix * Transpose(const BlockMatrix & A);
BlockMatrix * Mult(const BlockMatrix & A, const BlockMatrix & B);

inline void BlockMatrix::findGlobalRow(int iglobal, int & iblock, int & iloc) const
{
   if(iglobal > row_offsets[nRowBlocks])
      mfem_error("BlockMatrix::findGlobalRow");

   for(iblock = 0; iblock < nRowBlocks; ++iblock)
      if(row_offsets[iblock+1] > iglobal)
         break;

   iloc = iglobal - row_offsets[iblock];
}

inline void BlockMatrix::findGlobalCol(int jglobal, int & jblock, int & jloc) const
{
   if(jglobal > col_offsets[nColBlocks])
      mfem_error("BlockMatrix::findGlobalCol");

   for(jblock = 0; jblock < nColBlocks; ++jblock)
      if(col_offsets[jblock+1] > jglobal)
         break;

   jloc = jglobal - col_offsets[jblock];
}

}

#endif /* MFEM_BLOCKMATRIX */