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_MultiVectorFactory.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"

#include "Xpetra_MapExtractor.hpp"

#include "Xpetra_Matrix.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 = Matrix<>::scalar_type,
            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)
    {
      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(CrsMatrixFactory::Build(getRangeMap(r), npr, pftype));

      // Default view
      CreateDefaultView();
    }

#ifdef HAVE_XPETRA_THYRA

    BlockedCrsMatrix(const Teuchos::RCP<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 Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > 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,LO,GO,Node>::toXpetra(const_op);
            subRangeMaps[r] = xop->getRangeMap();
            break;
          }
        }
      }
      Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
      rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps));

      // 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);
      domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps));


      //std::cout << *(rangemaps_->getFullMap()) << std::endl;
      //Teuchos::RCP<const Xpetra::EpetraMap> epmap2 = Teuchos::rcp_dynamic_cast<const Xpetra::EpetraMap>(rangemaps_->getFullMap());
      //std::cout << epmap2->getEpetra_Map() << std::endl;

      //Teuchos::RCP<Thyra::LinearOpBase<Scalar> > fullOp = Teuchos::rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(thyraOp);

      //Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xFullOp =
      //                            Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toXpetra(fullOp);
      //std::cout << *(xFullOp->getRangeMap()) << std::endl;

      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) {
          if(thyraOp->blockExists(r,c)) {
            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);
            blocks_.push_back(xop);
          } else {
            // add empty block
            blocks_.push_back(CrsMatrixFactory::Build(getRangeMap(r), 0, Xpetra::DynamicProfile));
          }
        }
      }
      // Default view
      CreateDefaultView();
    }

  private:

    Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
      // 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);
        }
      }

      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) {
      throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
    }


    void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
      throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
    }

    void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map>& newMap) {
      throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
    }


    void replaceGlobalValues(GlobalOrdinal globalRow,
                             const ArrayView<const GlobalOrdinal> &cols,
                             const ArrayView<const Scalar>        &vals) {
      throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
    }


    void replaceLocalValues(LocalOrdinal localRow,
                            const ArrayView<const LocalOrdinal> &cols,
                            const ArrayView<const Scalar>       &vals) {
      throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
    }

    virtual void setAllToScalar(const Scalar& alpha) {
      throw Xpetra::Exceptions::RuntimeError("setAllToScalar not supported by BlockedCrsMatrix");
    }

    void scale(const Scalar& alpha) {
      throw Xpetra::Exceptions::RuntimeError("scale not supported by BlockedCrsMatrix");
    }




    void resumeFill(const RCP< ParameterList >& params = null) {
      throw Xpetra::Exceptions::RuntimeError("resumeFill not supported for block matrices");
    }

    void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
      throw Xpetra::Exceptions::RuntimeError("fillComplete with arguments not supported for block matrices");
    }

    void fillComplete(const RCP<ParameterList>& params = null) {
      for (size_t r = 0; r < Rows(); ++r)
        for (size_t c = 0; c < Cols(); ++c) {
          Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);

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

      // 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());

    }



    global_size_t getGlobalNumRows() const {
      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 {
      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 {
      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 {
      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 {
      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 {
      throw Xpetra::Exceptions::RuntimeError("getNumEntriesInLocalRow not supported by BlockedCrsMatrix");
    }


    global_size_t getGlobalNumDiags() const {
      throw Xpetra::Exceptions::RuntimeError("getGlobalNumDiags() not supported by BlockedCrsMatrix");
    }


    size_t getNodeNumDiags() const {
      throw Xpetra::Exceptions::RuntimeError("getNodeNumDiags() not supported by BlockedCrsMatrix");
    }


    size_t getGlobalMaxNumRowEntries() const {
      throw Xpetra::Exceptions::RuntimeError("getGlobalMaxNumRowEntries() not supported by BlockedCrsMatrix");
    }


    size_t getNodeMaxNumRowEntries() const {
      throw Xpetra::Exceptions::RuntimeError("getNodeMaxNumRowEntries() not supported by BlockedCrsMatrix");
    }


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


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

    bool isFillComplete() const {
      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 {
      throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
    }


    void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
      throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
    }


    void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
      throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
    }


    void getLocalDiagCopy(Vector& diag) const {
      throw Xpetra::Exceptions::RuntimeError("getLocalDiagCopy not supported by BlockedCrsMatrix" );
    }

    virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
      throw Xpetra::Exceptions::RuntimeError("getFrobeniusNorm() not supported by BlockedCrsMatrix, yet");
    }









    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
    {
      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." );

      RCP<const MultiVector> refX = rcpFromRef(X);
      RCP<MultiVector>       tmpY = MultiVectorFactory::Build(Y.getMap(), Y.getNumVectors());

      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());
          RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors());

          for (size_t col = 0; col < Cols(); col++) {
            RCP<const MultiVector> Xblock = domainmaps_->ExtractVector(refX, col);
            RCP<CrsMatrix>         Ablock = getMatrix(row, col);

            if (Ablock.is_null())
              continue;

            Ablock->apply(*Xblock, *tmpYblock);

            Yblock->update(one, *tmpYblock, one);
          }
          rangemaps_->InsertVector(Yblock, row, tmpY);
        }

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

          for (size_t row = 0; row<Rows(); row++) {
            RCP<const MultiVector> Xblock = rangemaps_->ExtractVector(refX, row);
            RCP<CrsMatrix>         Ablock = getMatrix(row, col);

            if (Ablock.is_null())
              continue;

            Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);

            Yblock->update(one, *tmpYblock, one);
          }
          domainmaps_->InsertVector(Yblock, col, tmpY);
        }
      }

      Y.update(alpha, *tmpY, beta);
    }

    RCP<const Map > getDomainMap() const            { return domainmaps_->getFullMap(); }

    RCP<const Map > getDomainMap(size_t i) const    { return domainmaps_->getMap(i); }

    RCP<const Map > getRangeMap() const             { return rangemaps_->getFullMap(); }

    RCP<const Map > getRangeMap(size_t i) const     { return rangemaps_->getMap(i); }

    RCP<const MapExtractor> getRangeMapExtractor()  { return rangemaps_; }

    RCP<const MapExtractor> getDomainMapExtractor() { return domainmaps_; }


    //{@

    const Teuchos::RCP< const Map > getMap() const {
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
    }

    void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
    }

    void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
    }

    void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
      throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
    }

    void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
      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;

        out << "fullRowMap" << std::endl;
        fullrowmap_->describe(out,verbLevel);

        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) {
          out << "Block(" << r << "," << c << ")" << std::endl;
          getMatrix(r,c)->describe(out,verbLevel);
        }
    }

    RCP<const CrsGraph> getCrsGraph() const {
      throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
    }




    virtual size_t Rows() const                                       { return rangemaps_->NumMaps(); }

    virtual size_t Cols() const                                       { return domainmaps_->NumMaps(); }

    Teuchos::RCP<CrsMatrix> getMatrix(size_t r, size_t c) const       { return blocks_[r*Cols()+c]; }

    void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
      // 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<CrsMatrix> Merge() const {
      using Teuchos::RCP;
      using Teuchos::rcp_dynamic_cast;
      Scalar one = ScalarTraits<SC>::one();

      RCP<CrsMatrix> sparse = CrsMatrixFactory::Build(fullrowmap_, 33);

      for (size_t i = 0; i < blocks_.size(); i++)
        if (blocks_[i] != Teuchos::null)
          this->Add(*blocks_[i], one, *sparse, one);

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

      return sparse;
    }

#ifdef HAVE_XPETRA_KOKKOS_REFACTOR
    typedef typename CrsMatrix::local_matrix_type local_matrix_type;

    local_matrix_type getLocalMatrix () const {
      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(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 CrsMatrix& A, const Scalar scalarA, CrsMatrix& B, const Scalar scalarB) const {
      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;

      size_t maxNumEntries = A.getNodeMaxNumRowEntries();

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

      RCP<const Map> rowMap = A.getRowMap();
      RCP<const Map> colMap = A.getColMap();

      ArrayView<const GO> rowGIDs = A.getRowMap()->getNodeElementList();
      for (size_t i = 0; i < A.getNodeNumRows(); i++) {
        GO row = rowGIDs[i];
        A.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() {
      // 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
    Teuchos::RCP<Map>                     fullrowmap_;        // full matrix    row map
    Teuchos::RCP<Map>                     fullcolmap_;        // full matrix column map

    std::vector<Teuchos::RCP<CrsMatrix> > blocks_;            // row major matrix block storage
#ifdef HAVE_XPETRA_THYRA
    Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > thyraOp_; 
#endif
};

} //namespace Xpetra

#define XPETRA_BLOCKEDCRSMATRIX_SHORT
#endif /* XPETRA_BLOCKEDCRSMATRIX_HPP */