Xpetra_BlockedCrsMatrix.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//             Xpetra: A linear algebra interface package
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Andrey Prokopenko (aprokop@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#ifndef XPETRA_BLOCKEDCRSMATRIX_HPP
#define XPETRA_BLOCKEDCRSMATRIX_HPP

#include <Kokkos_DefaultNode.hpp>

#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_Hashtable.hpp>

#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"

#include "Xpetra_MapFactory.hpp"
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_BlockedMultiVector.hpp"
#include "Xpetra_MultiVectorFactory.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"

#include "Xpetra_MapExtractor.hpp"

#include "Xpetra_Matrix.hpp"
#include "Xpetra_MatrixFactory.hpp"
#include "Xpetra_CrsMatrixWrap.hpp"

#ifdef HAVE_XPETRA_THYRA
#include <Thyra_ProductVectorSpaceBase.hpp>
#include <Thyra_VectorSpaceBase.hpp>
#include <Thyra_LinearOpBase.hpp>
#include <Thyra_BlockedLinearOpBase.hpp>
#include <Thyra_PhysicallyBlockedLinearOpBase.hpp>
#include "Xpetra_ThyraUtils.hpp"
#endif


namespace Xpetra {

  typedef std::string viewLabel_t;

  template <class Scalar,
            class LocalOrdinal,  // = typename Matrix<Scalar>::local_ordinal_type,
            class GlobalOrdinal, // = typename Matrix<Scalar, LocalOrdinal>::global_ordinal_type,
            class Node>          // = typename Matrix<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
  class BlockedCrsMatrix :
    public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
  public:
    typedef Scalar scalar_type;
    typedef LocalOrdinal local_ordinal_type;
    typedef GlobalOrdinal global_ordinal_type;
    typedef Node node_type;

  private:
#undef XPETRA_BLOCKEDCRSMATRIX_SHORT
#include "Xpetra_UseShortNames.hpp"

  public:




    BlockedCrsMatrix(Teuchos::RCP<const MapExtractor>& rangeMaps,
                     Teuchos::RCP<const MapExtractor>& domainMaps,
                     size_t npr, Xpetra::ProfileType pftype = Xpetra::DynamicProfile)
    : domainmaps_(domainMaps), rangemaps_(rangeMaps)
    {
      bRangeThyraMode_  = rangeMaps->getThyraMode();
      bDomainThyraMode_ = domainMaps->getThyraMode();

      blocks_.reserve(Rows()*Cols());

      // add CrsMatrix objects in row,column order
      for (size_t r = 0; r < Rows(); ++r)
        for (size_t c = 0; c < Cols(); ++c)
          blocks_.push_back(MatrixFactory::Build(getRangeMap(r,bRangeThyraMode_), npr, pftype));

      // Default view
      CreateDefaultView();
    }

#ifdef HAVE_XPETRA_THYRA

    BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
    : thyraOp_(thyraOp)
    {
      // extract information from Thyra blocked operator and rebuilt information
      const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace  = thyraOp->productRange();
      const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();

      int numRangeBlocks  = productRangeSpace->numBlocks();
      int numDomainBlocks = productDomainSpace->numBlocks();

      // build range map extractor from Thyra::BlockedLinearOpBase object
      std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
      for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
        for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
          if (thyraOp->blockExists(r,c)) {
            // we only need at least one block in each block row to extract the range map
            Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
            Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
                            Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(const_op);
            subRangeMaps[r] = xop->getRangeMap();
            break;
          }
        }
      }
      Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);

      // check whether the underlying Thyra operator uses Thyra-style numbering (default for most applications) or
      // Xpetra style numbering
      bool bRangeUseThyraStyleNumbering = false;
      size_t numAllElements = 0;
      for(size_t v = 0; v < subRangeMaps.size(); ++v) {
        numAllElements += subRangeMaps[v]->getGlobalNumElements();
      }
      if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering = true;
      rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps, bRangeUseThyraStyleNumbering));

      // build domain map extractor from Thyra::BlockedLinearOpBase object
      std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
      for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
        for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
          if (thyraOp->blockExists(r,c)) {
            // we only need at least one block in each block row to extract the range map
            Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
            Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
                            Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toXpetra(const_op);
            subDomainMaps[c] = xop->getDomainMap();
            break;
          }
        }
      }
      Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
      // plausibility check for numbering style (Xpetra or Thyra)
      bool bDomainUseThyraStyleNumbering = false;
      numAllElements = 0;
      for(size_t v = 0; v < subDomainMaps.size(); ++v) {
        numAllElements += subDomainMaps[v]->getGlobalNumElements();
      }
      if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering = true;
      domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps, bDomainUseThyraStyleNumbering));

      // store numbering mode
      bRangeThyraMode_  = bRangeUseThyraStyleNumbering;
      bDomainThyraMode_ = bDomainUseThyraStyleNumbering;

      // add CrsMatrix objects in row,column order
      blocks_.reserve(Rows()*Cols());
      for (size_t r = 0; r < Rows(); ++r) {
        for (size_t c = 0; c < Cols(); ++c) {
          if(thyraOp->blockExists(r,c)) {
            // TODO we do not support nested Thyra operators here!
            Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
            Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
            Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
                Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toXpetra(op);
            Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node> > xwrap =
                Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node>(xop));
            blocks_.push_back(xwrap);
          } else {
            // add empty block
            blocks_.push_back(MatrixFactory::Build(getRangeMap(r,bRangeThyraMode_), 0, Xpetra::DynamicProfile));
          }
        }
      }
      // Default view
      CreateDefaultView();
    }

  private:

    Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
      // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor

      // merge submaps to global map
      std::vector<GlobalOrdinal> gids;
      for(size_t tt = 0; tt<subMaps.size(); ++tt) {
        Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > subMap = subMaps[tt];
        for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(subMap->getNodeNumElements()); ++l) {
          GlobalOrdinal gid = subMap->getGlobalElement(l);
          gids.push_back(gid);
        }
      }

      // we have to sort the matrix entries and get rid of the double entries
      // since we use this to detect Thyra-style numbering or Xpetra-style
      // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
      // the correct row maps.
      const GO INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
      std::sort(gids.begin(), gids.end());
      gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
      Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
      Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(subMaps[0]->lib(), INVALID, gidsView, subMaps[0]->getIndexBase(), subMaps[0]->getComm());
      return fullMap;
    }

  public:
#endif

    virtual ~BlockedCrsMatrix() {}






    void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertGlobalValues");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
    }


    void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertLocalValues");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
    }

    void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map>& newMap) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
    }


    void replaceGlobalValues(GlobalOrdinal globalRow,
                             const ArrayView<const GlobalOrdinal> &cols,
                             const ArrayView<const Scalar>        &vals) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceGlobalValues");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
    }


    void replaceLocalValues(LocalOrdinal localRow,
                            const ArrayView<const LocalOrdinal> &cols,
                            const ArrayView<const Scalar>       &vals) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceLocalValues");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
    }

    virtual void setAllToScalar(const Scalar& alpha) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::setAllToScalar");
      for (size_t row = 0; row < Rows(); row++) {
        for (size_t col = 0; col < Cols(); col++) {
          if (!getMatrix(row,col).is_null()) {
            getMatrix(row,col)->setAllToScalar(alpha);
          }
        }
      }
    }

    void scale(const Scalar& alpha) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::scale");
      for (size_t row = 0; row < Rows(); row++) {
        for (size_t col = 0; col < Cols(); col++) {
          if (!getMatrix(row,col).is_null()) {
            getMatrix(row,col)->scale(alpha);
          }
        }
      }
    }




    void resumeFill(const RCP< ParameterList >& params = null) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::resumeFill");
      for (size_t row = 0; row < Rows(); row++) {
        for (size_t col = 0; col < Cols(); col++) {
          if (!getMatrix(row,col).is_null()) {
            getMatrix(row,col)->resumeFill(params);
          }
        }
      }
    }

    void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
        return;
      }
      fillComplete(params);
    }

    void fillComplete(const RCP<ParameterList>& params = null) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
      TEUCHOS_TEST_FOR_EXCEPTION(rangemaps_==Teuchos::null, Xpetra::Exceptions::RuntimeError,"BlockedCrsMatrix::fillComplete: rangemaps_ is not set. Error.");

      for (size_t r = 0; r < Rows(); ++r)
        for (size_t c = 0; c < Cols(); ++c) {
          if(getMatrix(r,c) != Teuchos::null) {
            Teuchos::RCP<Matrix> Ablock = getMatrix(r,c);

            if (Ablock != Teuchos::null && !Ablock->isFillComplete())
              Ablock->fillComplete(getDomainMap(c, bDomainThyraMode_), getRangeMap(r, bRangeThyraMode_), params);
          }
        }

#if 0
      // get full row map
      RCP<const Map> rangeMap = rangemaps_->getFullMap();
      fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());

      // build full col map
      fullcolmap_ = Teuchos::null; // delete old full column map

      std::vector<GO> colmapentries;
      for (size_t c = 0; c < Cols(); ++c) {
        // copy all local column lids of all block rows to colset
        std::set<GO> colset;
        for (size_t r = 0; r < Rows(); ++r) {
          Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);

          if (Ablock != Teuchos::null) {
            Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
            Teuchos::RCP<const Map>      colmap      = Ablock->getColMap();
            copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
          }
        }

        // remove duplicates (entries which are in column maps of more than one block row)
        colmapentries.reserve(colmapentries.size() + colset.size());
        copy(colset.begin(), colset.end(), back_inserter(colmapentries));
        sort(colmapentries.begin(), colmapentries.end());
        typename std::vector<GO>::iterator gendLocation;
        gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
        colmapentries.erase(gendLocation,colmapentries.end());
      }

      // sum up number of local elements
      size_t numGlobalElements = 0;
      Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));

      // store global full column map
      const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
      fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
#endif
    }



    global_size_t getGlobalNumRows() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumRows");
      global_size_t globalNumRows = 0;

      for (size_t row = 0; row < Rows(); row++)
        for (size_t col = 0; col < Cols(); col++)
          if (!getMatrix(row,col).is_null()) {
            globalNumRows += getMatrix(row,col)->getGlobalNumRows();
            break; // we need only one non-null matrix in a row
          }

      return globalNumRows;
    }


    global_size_t getGlobalNumCols() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumCols");
      global_size_t globalNumCols = 0;

      for (size_t col = 0; col < Cols(); col++)
        for (size_t row = 0; row < Rows(); row++)
          if (!getMatrix(row,col).is_null()) {
            globalNumCols += getMatrix(row,col)->getGlobalNumCols();
            break; // we need only one non-null matrix in a col
          }

      return globalNumCols;
    }

    size_t getNodeNumRows() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumRows");
      global_size_t nodeNumRows = 0;

      for (size_t row = 0; row < Rows(); ++row)
        for (size_t col = 0; col < Cols(); col++)
          if (!getMatrix(row,col).is_null()) {
            nodeNumRows += getMatrix(row,col)->getNodeNumRows();
            break; // we need only one non-null matrix in a row
          }

      return nodeNumRows;
    }

    global_size_t getGlobalNumEntries() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumEntries");
      global_size_t globalNumEntries = 0;

      for (size_t row = 0; row < Rows(); ++row)
        for (size_t col = 0; col < Cols(); ++col)
          if (!getMatrix(row,col).is_null())
            globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();

      return globalNumEntries;
    }

    size_t getNodeNumEntries() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumEntries");
      global_size_t nodeNumEntries = 0;

      for (size_t row = 0; row < Rows(); ++row)
        for (size_t col = 0; col < Cols(); ++col)
          if (!getMatrix(row,col).is_null())
            nodeNumEntries += getMatrix(row,col)->getNodeNumEntries();

      return nodeNumEntries;
    }


    size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getNumEntriesInLocalRow(localRow);
      }
      throw Xpetra::Exceptions::RuntimeError("getNumEntriesInLocalRow not supported by BlockedCrsMatrix");
    }


    global_size_t getGlobalNumDiags() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumDiags");
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getGlobalNumDiags();
      }
      throw Xpetra::Exceptions::RuntimeError("getGlobalNumDiags() not supported by BlockedCrsMatrix");
    }


    size_t getNodeNumDiags() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumDiags");
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getNodeNumDiags();
      }
      throw Xpetra::Exceptions::RuntimeError("getNodeNumDiags() not supported by BlockedCrsMatrix");
    }


    size_t getGlobalMaxNumRowEntries() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
      global_size_t globalMaxEntries = 0;

      for (size_t row = 0; row < Rows(); row++) {
        global_size_t globalMaxEntriesBlockRows = 0;
        for (size_t col = 0; col < Cols(); col++) {
          if (!getMatrix(row,col).is_null()) {
            globalMaxEntriesBlockRows += getMatrix(row,col)->getGlobalMaxNumRowEntries();
          }
        }
        if(globalMaxEntriesBlockRows > globalMaxEntries)
          globalMaxEntries = globalMaxEntriesBlockRows;
      }
      return globalMaxEntries;
    }


    size_t getNodeMaxNumRowEntries() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeMaxNumRowEntries");
      size_t localMaxEntries = 0;

      for (size_t row = 0; row < Rows(); row++) {
        size_t localMaxEntriesBlockRows = 0;
        for (size_t col = 0; col < Cols(); col++) {
          if (!getMatrix(row,col).is_null()) {
            localMaxEntriesBlockRows += getMatrix(row,col)->getNodeMaxNumRowEntries();
          }
        }
        if(localMaxEntriesBlockRows > localMaxEntries)
          localMaxEntries = localMaxEntriesBlockRows;
      }
      return localMaxEntries;
    }


    bool isLocallyIndexed() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::isLocallyIndexed");
      for (size_t i = 0; i < blocks_.size(); ++i)
        if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
          return false;
      return true;
    }


    bool isGloballyIndexed() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::isGloballyIndexed");
      for (size_t i = 0; i < blocks_.size(); i++)
        if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
          return false;
      return true;
    }

    bool isFillComplete() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::isFillComplete");
      for (size_t i = 0; i < blocks_.size(); i++)
        if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
          return false;
      return true;
    }


    virtual void getLocalRowCopy(LocalOrdinal LocalRow,
                                 const ArrayView<LocalOrdinal>& Indices,
                                 const ArrayView<Scalar>& Values,
                                 size_t &NumEntries) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowCopy");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
    }


    void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalRowView");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
    }


    void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowView");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
    }


    void getLocalDiagCopy(Vector& diag) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalDiagCopy");
      XPETRA_TEST_FOR_EXCEPTION(diag.getMap()->isSameAs(*rangemaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
        "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the full map of the blocked operator." );

      TEUCHOS_TEST_FOR_EXCEPTION(Rows() != Cols(), Xpetra::Exceptions::RuntimeError,
        "BlockedCrsMatrix::getLocalDiagCopy(): you cannot extract the diagonal of a "<< Rows() << "x"<<Cols()<<" blocked matrix." );

      for (size_t row = 0; row < Rows(); ++row) {
        if (!getMatrix(row,row).is_null()) {
          // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
          bool bThyraMode = rangemaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,row)) == Teuchos::null);
          RCP<Vector> dd = rangemaps_->getVector(row, bThyraMode);
          getMatrix(row,row)->getLocalDiagCopy(*dd);
          rangemaps_->InsertVector(*dd,row,diag,bThyraMode);
        }
      }
    }

    void leftScale (const Vector& x) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::leftScale");
      XPETRA_TEST_FOR_EXCEPTION(x.getMap()->isSameAs(*rangemaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
        "BlockedCrsMatrix::leftScale(): the map of the vector x is not compatible with the full map of the blocked operator." );

      RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);

      for (size_t row = 0; row < Rows(); ++row) {
        for (size_t col = 0; col < Cols(); ++col) {
          if(getMatrix(row,col)!=Teuchos::null) {
            // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
            bool bThyraMode = rangemaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,col)) == Teuchos::null);
            RCP<Vector> xx = rangemaps_->ExtractVector(rcpx,row,bThyraMode);
            getMatrix(row,col)->leftScale(*xx);
          }
        }
      }
    }

    void rightScale (const Vector& x) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::rightScale");
      XPETRA_TEST_FOR_EXCEPTION(x.getMap()->isSameAs(*domainmaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
        "BlockedCrsMatrix::rightScale(): the map of the vector x is not compatible with the full map of the blocked operator." );

      RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);

      for (size_t col = 0; col < Cols(); ++col) {
        for (size_t row = 0; row < Rows(); ++row) {
          if(getMatrix(row,col)!=Teuchos::null) {
            // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
            bool bThyraMode = domainmaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,col)) == Teuchos::null);
            RCP<Vector> xx = domainmaps_->ExtractVector(rcpx,col,bThyraMode);
            getMatrix(row,col)->rightScale(*xx);
          }
        }
      }
    }


    virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFrobeniusNorm");
      typename ScalarTraits<Scalar>::magnitudeType ret = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero());
      for (size_t col = 0; col < Cols(); ++col) {
        for (size_t row = 0; row < Rows(); ++row) {
          if(getMatrix(row,col)!=Teuchos::null) {
            typename ScalarTraits<Scalar>::magnitudeType n = getMatrix(row,col)->getFrobeniusNorm();
            ret += n * n;
          }
        }
      }
      return Teuchos::ScalarTraits< typename ScalarTraits<Scalar>::magnitudeType >::squareroot(ret);
    }









    virtual void apply(const MultiVector& X, MultiVector& Y,
                       Teuchos::ETransp mode = Teuchos::NO_TRANS,
                       Scalar alpha = ScalarTraits<Scalar>::one(),
                       Scalar beta  = ScalarTraits<Scalar>::zero()) const
    {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::apply");
      //using Teuchos::RCP;

      TEUCHOS_TEST_FOR_EXCEPTION(mode != Teuchos::NO_TRANS && mode != Teuchos::TRANS, Xpetra::Exceptions::RuntimeError,
                                 "apply() only supports the following modes: NO_TRANS and TRANS." );

      // check whether input parameters are blocked or not
      RCP<const MultiVector>         refX = rcpFromRef(X);
      RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
      RCP<MultiVector>               tmpY = rcpFromRef(Y);
      RCP<BlockedMultiVector>       tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);

      bool bBlockedX = (refbX != Teuchos::null) ? true : false;
      bool bBlockedY = (tmpbY != Teuchos::null) ? true : false;

      // create (temporary) vectors for output
      // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
      tmpY = MultiVectorFactory::Build(Y.getMap(), Y.getNumVectors(), true);
      if (bBlockedY == true) tmpbY = Teuchos::rcp(new BlockedMultiVector(rangemaps_,tmpY));

      SC one = ScalarTraits<SC>::one();

      if (mode == Teuchos::NO_TRANS) {

        for (size_t row = 0; row < Rows(); row++) {
          RCP<MultiVector>    Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
          for (size_t col = 0; col < Cols(); col++) {

            // extract matrix block
            RCP<Matrix>            Ablock = getMatrix(row, col);

            if (Ablock.is_null())
              continue;

            // check whether Ablock is itself a blocked operator
            // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
            bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;

            // input/output vectors for local block operation
            RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A

            // extract sub part of X using Xpetra or Thyra GIDs
            // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
            // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
            if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
            else          Xblock = domainmaps_->ExtractVector(refX,  col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
            RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bBlockedSubMatrix == true ? false : bRangeThyraMode_, false);  // subpart of Y vector containing part of solution of Xblock applied to Ablock
            Ablock->apply(*Xblock, *tmpYblock);

            // If Ablock is a blocked operator the local vectors are using (pseudo) Xpetra-style gids
            // that have to be translated to Thyra based GIDs if bRangeThyraMode is set
            if(bBlockedSubMatrix == true && bRangeThyraMode_ == true) {
              tmpYblock->replaceMap(rangemaps_->getMap(row, true)); // switch to Thyra maps (compatible to Yblock)
            }
            Yblock->update(one, *tmpYblock, one);
          }
          if(bBlockedY) rangemaps_->InsertVector(Yblock, row, tmpbY, bRangeThyraMode_);
          else          rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
        }

      } else if (mode == Teuchos::TRANS) {
        // TODO: test me!
        for (size_t col = 0; col < Cols(); col++) {
          RCP<MultiVector>    Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, true);

          for (size_t row = 0; row<Rows(); row++) {
            RCP<Matrix>            Ablock = getMatrix(row, col);

            if (Ablock.is_null())
              continue;

            // check whether Ablock is itself a blocked operator
            // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
            bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;

            RCP<const MultiVector> Xblock = Teuchos::null;

            // extract sub part of X using Xpetra or Thyra GIDs
            if(bBlockedX) Xblock = rangemaps_->ExtractVector(refbX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
            else          Xblock = rangemaps_->ExtractVector(refX,  row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
            RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bBlockedSubMatrix == true ? false : bDomainThyraMode_, false);
            Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);

            // If Ablock is a blocked operator the local vectors are using (pseudo) Xpetra-style gids
            // that have to be translated to Thyra based GIDs if bRangeThyraMode is set
            if(bBlockedSubMatrix == true && bDomainThyraMode_ == true) {
              tmpYblock->replaceMap(domainmaps_->getMap(col, true)); // switch to Thyra maps (compatible to Yblock)
            }

            Yblock->update(one, *tmpYblock, one);
          }
          if(bBlockedY) domainmaps_->InsertVector(Yblock, col, tmpbY, bDomainThyraMode_);
          else          domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
        }
      }
      if(bBlockedY) Y.update(alpha, *tmpbY, beta);
      else          Y.update(alpha, *tmpY, beta);
    }

    RCP<const Map > getDomainMap() const            { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap()"); return domainmaps_->getFullMap(); }

    RCP<const Map > getDomainMap(size_t i) const    { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t)"); return domainmaps_->getMap(i, bDomainThyraMode_); }

    RCP<const Map > getDomainMap(size_t i, bool bThyraMode) const    { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t,bool)"); return domainmaps_->getMap(i, bThyraMode); }

    RCP<const Map > getRangeMap() const             { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getFullMap(); }

    RCP<const Map > getRangeMap(size_t i) const     { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t)"); return rangemaps_->getMap(i, bRangeThyraMode_); }

    RCP<const Map > getRangeMap(size_t i, bool bThyraMode) const     { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t,bool)"); return rangemaps_->getMap(i, bThyraMode); }

    RCP<const MapExtractor> getRangeMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMapExtractor()"); return rangemaps_; }

    RCP<const MapExtractor> getDomainMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMapExtractor()"); return domainmaps_; }


    //{@

    const Teuchos::RCP< const Map > getMap() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMap");
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getMap();
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
    }

    void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->doImport(source, importer, CM);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
    }

    void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->doExport(dest, importer, CM);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
    }

    void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->doImport(source, exporter, CM);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
    }

    void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
      if (Rows() == 1 && Cols () == 1) {
        getMatrix(0,0)->doExport(dest, exporter, CM);
        return;
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
    }

    // @}



    std::string description() const                       { return "Xpetra_BlockedCrsMatrix.description()"; }

    void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
      out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;

      if (isFillComplete()) {
        out << "BlockMatrix is fillComplete" << std::endl;

        /*if(fullrowmap_ != Teuchos::null) {
          out << "fullRowMap" << std::endl;
          fullrowmap_->describe(out,verbLevel);
        } else {
          out << "fullRowMap not set. Check whether block matrix is properly fillCompleted!" << std::endl;
        }*/

        //out << "fullColMap" << std::endl;
        //fullcolmap_->describe(out,verbLevel);

      } else {
        out << "BlockMatrix is NOT fillComplete" << std::endl;
      }

      for (size_t r = 0; r < Rows(); ++r)
        for (size_t c = 0; c < Cols(); ++c) {
          if(getMatrix(r,c)!=Teuchos::null) {
            out << "Block(" << r << "," << c << ")" << std::endl;
            getMatrix(r,c)->describe(out,verbLevel);
          } else out << "Block(" << r << "," << c << ") = null" << std::endl;
        }
    }

    RCP<const CrsGraph> getCrsGraph() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsGraph");
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getCrsGraph();
      }
      throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
    }




    virtual size_t Rows() const                                       { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Rows"); return rangemaps_->NumMaps(); }

    virtual size_t Cols() const                                       { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Cols"); return domainmaps_->NumMaps(); }

    Teuchos::RCP<Matrix> getCrsMatrix() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsMatrix");
      TEUCHOS_TEST_FOR_EXCEPTION(Rows()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Rows() << " block rows, though.");
      TEUCHOS_TEST_FOR_EXCEPTION(Cols()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Cols() << " block columns, though.");

      RCP<Matrix> mat = getMatrix(0,0);
      RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
      if (bmat == Teuchos::null) return mat;
      return bmat->getCrsMatrix();
    }

    Teuchos::RCP<Matrix> getInnermostCrsMatrix() {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getInnermostCrsMatrix");
      size_t row = Rows()+1, col = Cols()+1;
      for (size_t r = 0; r < Rows(); ++r)
        for(size_t c = 0; c < Cols(); ++c)
          if (getMatrix(r,c) != Teuchos::null) {
            row = r;
            col = c;
            break;
          }
      TEUCHOS_TEST_FOR_EXCEPTION(row == Rows()+1 || col == Cols()+1, Xpetra::Exceptions::Incompatible, "Xpetra::BlockedCrsMatrix::getInnermostCrsMatrix: Could not find a non-zero sub-block in blocked operator.")
      RCP<Matrix> mm = getMatrix(row,col);
      RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mm);
      if (bmat == Teuchos::null) return mm;
      return bmat->getInnermostCrsMatrix();
    }

    Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const       {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMatrix");
      TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
      TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");

      // transfer strided/blocked map information
      if (blocks_[r*Cols()+c] != Teuchos::null &&
          blocks_[r*Cols()+c]->IsView("stridedMaps") == false)
        blocks_[r*Cols()+c]->CreateView("stridedMaps", getRangeMap(r,bRangeThyraMode_), getDomainMap(c,bDomainThyraMode_));
      return blocks_[r*Cols()+c];
    }

    //void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
    void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat) {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::setMatrix");
      // TODO: if filled -> return error

      TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
      TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");

      // set matrix
      blocks_[r*Cols() + c] = mat;
    }

    // NOTE: This is a rather expensive operation, since all blocks are copied
    // into a new big CrsMatrix
    Teuchos::RCP<Matrix> Merge() const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::Merge");
      using Teuchos::RCP;
      using Teuchos::rcp_dynamic_cast;
      Scalar one = ScalarTraits<SC>::one();

      TEUCHOS_TEST_FOR_EXCEPTION(bRangeThyraMode_ != bDomainThyraMode_, Xpetra::Exceptions::RuntimeError,
                                 "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );

      TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, Xpetra::Exceptions::RuntimeError,
                                 "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );

      RCP<Matrix> sparse = MatrixFactory::Build(getRangeMapExtractor()->getFullMap(), 33);

      if(bRangeThyraMode_ == false) {
        // Xpetra mode
        for (size_t i = 0; i < Rows(); i++) {
          for (size_t j = 0; j < Cols(); j++) {
            if (getMatrix(i,j) != Teuchos::null) {
              RCP<const Matrix> mat = getMatrix(i,j);

              // recursively call Merge routine
              RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
              if (bMat != Teuchos::null) mat = bMat->Merge();

              bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
              TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
                                         "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );

              // jump over empty blocks
              if(mat->getNodeNumEntries() == 0) continue;

              this->Add(*mat, one, *sparse, one);
            }
          }
        }
      } else {
        // Thyra mode
        for (size_t i = 0; i < Rows(); i++) {
          for (size_t j = 0; j < Cols(); j++) {
            if (getMatrix(i,j) != Teuchos::null) {
              RCP<const Matrix> mat = getMatrix(i,j);
              // recursively call Merge routine
              RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
              if (bMat != Teuchos::null) mat = bMat->Merge();

              bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
              TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
                                         "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );

              // check whether we have a CrsMatrix block (no blocked operator)
              RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(mat);
              TEUCHOS_ASSERT(crsMat != Teuchos::null);

              // these are the thyra style maps of the matrix
              RCP<const Map> trowMap = mat->getRowMap();
              RCP<const Map> tcolMap = mat->getColMap();
              RCP<const Map> tdomMap = mat->getDomainMap();

              // get Xpetra maps
              RCP<const Map> xrowMap = getRangeMapExtractor()->getMap(i,false);
              RCP<const Map> xdomMap = getDomainMapExtractor()->getMap(j,false);

              // generate column map with Xpetra GIDs
              // We have to do this separately for each block since the column
              // map of each block might be different in the same block column
              Teuchos::RCP<Map> xcolMap = MapUtils::transformThyra2XpetraGIDs(
                    *tcolMap,
                    *tdomMap,
                    *xdomMap);

              // jump over empty blocks
              if(mat->getNodeNumEntries() == 0) continue;

              size_t maxNumEntries = mat->getNodeMaxNumRowEntries();

              size_t    numEntries;
              Array<GO> inds (maxNumEntries);
              Array<GO> inds2(maxNumEntries);
              Array<SC> vals (maxNumEntries);

              // loop over all rows and add entries
              for (size_t k = 0; k < mat->getNodeNumRows(); k++) {
                GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
                crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);

                // create new indices array
                for (size_t l = 0; l < numEntries; ++l) {
                  LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
                  inds2[l] = xcolMap->getGlobalElement(lid);
                }

                GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
                sparse->insertGlobalValues(
                    rowXGID, inds2(0, numEntries),
                    vals(0, numEntries));
              }
            }
          }
        }
      }

      sparse->fillComplete(getDomainMap(), getRangeMap());

      TEUCHOS_TEST_FOR_EXCEPTION(sparse->getNodeNumEntries() != getNodeNumEntries(), Xpetra::Exceptions::RuntimeError,
                                 "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );

      TEUCHOS_TEST_FOR_EXCEPTION(sparse->getGlobalNumEntries() != getGlobalNumEntries(), Xpetra::Exceptions::RuntimeError,
                                 "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );

      return sparse;
    }

#ifdef HAVE_XPETRA_KOKKOS_REFACTOR
    typedef typename CrsMatrix::local_matrix_type local_matrix_type;

    local_matrix_type getLocalMatrix () const {
      if (Rows() == 1 && Cols () == 1) {
        return getMatrix(0,0)->getLocalMatrix();
      }
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
    }
#endif

#ifdef HAVE_XPETRA_THYRA
    Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > getThyraOperator() {
      Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thOp =
          Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*this));
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thOp));

      Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > thbOp =
          Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thbOp));
      return thbOp;
    }
#endif

  private:



    void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::Add");
      TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError,
                                 "Matrix A is not completed");
      using Teuchos::Array;
      using Teuchos::ArrayView;

      B.scale(scalarB);

      Scalar one  = ScalarTraits<SC>::one();
      Scalar zero = ScalarTraits<SC>::zero();

      if (scalarA == zero)
        return;

      Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
      Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(rcpA);
      TEUCHOS_TEST_FOR_EXCEPTION(rcpAwrap == Teuchos::null, Xpetra::Exceptions::BadCast,
                                       "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
      Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();

      size_t maxNumEntries = crsA->getNodeMaxNumRowEntries();

      size_t    numEntries;
      Array<GO> inds(maxNumEntries);
      Array<SC> vals(maxNumEntries);

      RCP<const Map> rowMap = crsA->getRowMap();
      RCP<const Map> colMap = crsA->getColMap();

      ArrayView<const GO> rowGIDs = crsA->getRowMap()->getNodeElementList();
      for (size_t i = 0; i < crsA->getNodeNumRows(); i++) {
        GO row = rowGIDs[i];
        crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);

        if (scalarA != one)
          for (size_t j = 0; j < numEntries; ++j)
            vals[j] *= scalarA;

        B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
      }
    }


    // Default view is created after fillComplete()
    // Because ColMap might not be available before fillComplete().
    void CreateDefaultView() {
      XPETRA_MONITOR("XpetraBlockedCrsMatrix::CreateDefaultView");

      // Create default view
      this->defaultViewLabel_ = "point";
      this->CreateView(this->GetDefaultViewLabel(), getRangeMap(), getDomainMap());

      // Set current view
      this->currentViewLabel_ = this->GetDefaultViewLabel();
    }

  private:
    Teuchos::RCP<const MapExtractor>      domainmaps_;        // full        domain map together with all partial domain maps
    Teuchos::RCP<const MapExtractor>      rangemaps_;         // full         range map together with all partial domain maps

    std::vector<Teuchos::RCP<Matrix> > blocks_;            // row major matrix block storage
#ifdef HAVE_XPETRA_THYRA
    Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > thyraOp_; 
#endif
    bool                                   bRangeThyraMode_;  
    bool                                   bDomainThyraMode_; 

};

} //namespace Xpetra

#define XPETRA_BLOCKEDCRSMATRIX_SHORT
#endif /* XPETRA_BLOCKEDCRSMATRIX_HPP */