Xpetra_Matrix.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
 
// WARNING: This code is experimental. Backwards compatibility should not be expected.
 
#ifndef XPETRA_MATRIX_HPP
#define XPETRA_MATRIX_HPP
 
#include <Kokkos_DefaultNode.hpp>
 
#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"
 
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"
#include "Xpetra_MatrixView.hpp"
#include "Xpetra_Operator.hpp"
#include "Xpetra_StridedMap.hpp"
#include "Xpetra_StridedMapFactory.hpp"
 
#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_Hashtable.hpp>
 
namespace Xpetra {
 
  typedef std::string viewLabel_t;
 
  template <class Scalar,
            class LocalOrdinal,
            class GlobalOrdinal,
            class Node = KokkosClassic::DefaultNode::DefaultNodeType>
  class Matrix : public Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > {
    typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> Map;
    typedef Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> CrsMatrix;
    typedef Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node> CrsGraph;
    typedef Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node> CrsMatrixFactory;
    typedef Xpetra::MatrixView<Scalar, LocalOrdinal, GlobalOrdinal, Node> MatrixView;
 
  public:
    typedef Scalar          scalar_type;
    typedef LocalOrdinal    local_ordinal_type;
    typedef GlobalOrdinal   global_ordinal_type;
    typedef Node            node_type;
 
#ifdef HAVE_XPETRA_KOKKOS_REFACTOR
#ifdef HAVE_XPETRA_TPETRA
    typedef typename CrsMatrix::local_matrix_type local_matrix_type;
#endif
#endif
 
 
 
    Matrix() { }
 
    virtual ~Matrix() { }
 
 
 
    void CreateView(viewLabel_t viewLabel, const RCP<const Map> & rowMap, const RCP<const Map> & colMap) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == true, Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.CreateView(): a view labeled '" + viewLabel + "' already exist.");
      RCP<MatrixView> view = rcp(new MatrixView(rowMap, colMap));
      operatorViewTable_.put(viewLabel, view);
    }
 
    // JG TODO: why this is a member function??
    void CreateView(const viewLabel_t viewLabel, const RCP<const Matrix>& A, bool transposeA = false, const RCP<const Matrix>& B = Teuchos::null, bool transposeB = false) {
      RCP<const Map> domainMap = Teuchos::null;
      RCP<const Map> rangeMap  = Teuchos::null;
 
      const size_t        blkSize = 1;
      std::vector<size_t> stridingInfo(1, blkSize);
      LocalOrdinal        stridedBlockId = -1;
 
 
      if (A->IsView(viewLabel)) {
        rangeMap  = transposeA ? A->getColMap(viewLabel) : A->getRowMap(viewLabel);
        domainMap = transposeA ? A->getRowMap(viewLabel) : A->getColMap(viewLabel); // will be overwritten if B != Teuchos::null
 
      } else {
        rangeMap  = transposeA ? A->getDomainMap()       : A->getRangeMap();
        domainMap = transposeA ? A->getRangeMap()        : A->getDomainMap();
 
        if (viewLabel == "stridedMaps") {
          rangeMap  = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(rangeMap,  stridingInfo, stridedBlockId);
          domainMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(domainMap, stridingInfo, stridedBlockId);
        }
      }
 
      if (B != Teuchos::null ) {
        // B has strided Maps
 
        if (B->IsView(viewLabel)) {
          domainMap = transposeB ? B->getRowMap(viewLabel) : B->getColMap(viewLabel);
 
        } else {
          domainMap = transposeB ? B->getRangeMap()        : B->getDomainMap();
 
          if (viewLabel == "stridedMaps")
            domainMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(domainMap, stridingInfo, stridedBlockId);
        }
      }
 
 
      if (IsView(viewLabel))
        RemoveView(viewLabel);
 
      CreateView(viewLabel, rangeMap, domainMap);
    }
 
    void PrintViews(Teuchos::FancyOStream &out) const {
      int last = out.getOutputToRootOnly();
      Teuchos::OSTab tab(out);
      out.setOutputToRootOnly(0);
      Teuchos::Array<viewLabel_t> viewLabels;
      Teuchos::Array<RCP<MatrixView> > viewList;
      operatorViewTable_.arrayify(viewLabels,viewList);
      out << "views associated with this operator" << std::endl;
      for (int i=0; i<viewLabels.size(); ++i)
        out << viewLabels[i] << std::endl;
      out.setOutputToRootOnly(last);
    }
 
 
    void RemoveView(const viewLabel_t viewLabel) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.RemoveView(): view '" + viewLabel + "' does not exist.");
      TEUCHOS_TEST_FOR_EXCEPTION(viewLabel == GetDefaultViewLabel(), Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.RemoveView(): view '" + viewLabel + "' is the default view and cannot be removed.");
      operatorViewTable_.remove(viewLabel);
    }
 
    const viewLabel_t SwitchToView(const viewLabel_t viewLabel) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.SwitchToView(): view '" + viewLabel + "' does not exist.");
      viewLabel_t oldViewLabel = GetCurrentViewLabel();
      currentViewLabel_ = viewLabel;
      return oldViewLabel;
    }
 
    bool IsView(const viewLabel_t viewLabel) const {
      return operatorViewTable_.containsKey(viewLabel);
    }
 
    const viewLabel_t SwitchToDefaultView() { return SwitchToView(GetDefaultViewLabel()); }
 
    const viewLabel_t & GetDefaultViewLabel() const { return defaultViewLabel_; }
 
    const viewLabel_t & GetCurrentViewLabel() const { return currentViewLabel_; }
 
 
 
 
 
    virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals) = 0;
 
 
    virtual void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals) = 0;
 
 
    virtual void replaceGlobalValues(GlobalOrdinal globalRow,
                                     const ArrayView<const GlobalOrdinal> &cols,
                                     const ArrayView<const Scalar>        &vals) = 0;
 
 
    virtual void replaceLocalValues(LocalOrdinal localRow,
                                    const ArrayView<const LocalOrdinal> &cols,
                                    const ArrayView<const Scalar>       &vals) = 0;
 
    virtual void setAllToScalar(const Scalar &alpha)= 0;
 
    virtual void scale(const Scalar &alpha)= 0;
 
 
 
 
    virtual void resumeFill(const RCP< ParameterList > &params=null) = 0;
 
    virtual void fillComplete(const RCP<const Map> &domainMap, const RCP<const Map> &rangeMap, const RCP<ParameterList> &params = null) =0;
 
    //TODO : Get ride of "Tpetra"::OptimizeOption
    virtual void fillComplete(const RCP<ParameterList> &params=null) =0;
 
 
 
 
    virtual const RCP<const Map> & getRowMap() const { return getRowMap(GetCurrentViewLabel()); }
 
    virtual const RCP<const Map> & getRowMap(viewLabel_t viewLabel) const {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.GetRowMap(): view '" + viewLabel + "' does not exist.");
      return operatorViewTable_.get(viewLabel)->GetRowMap();
    }
 
    virtual const RCP<const Map> & getColMap() const { return getColMap(GetCurrentViewLabel()); }
 
    virtual const RCP<const Map> & getColMap(viewLabel_t viewLabel) const {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Matrix.GetColMap(): view '" + viewLabel + "' does not exist.");
      return operatorViewTable_.get(viewLabel)->GetColMap();
    }
 
 
    virtual global_size_t getGlobalNumRows() const =0;
 
 
    virtual global_size_t getGlobalNumCols() const =0;
 
    virtual size_t getNodeNumRows() const =0;
 
    virtual global_size_t getGlobalNumEntries() const =0;
 
    virtual size_t getNodeNumEntries() const =0;
 
 
    virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0;
 
 
    virtual size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const =0;
 
 
    virtual size_t getGlobalMaxNumRowEntries() const =0;
 
 
    virtual size_t getNodeMaxNumRowEntries() const =0;
 
    virtual bool isLocallyIndexed() const =0;
 
    virtual bool isGloballyIndexed() const =0;
 
    virtual bool isFillComplete() const =0;
 
 
    virtual void getLocalRowCopy(LocalOrdinal LocalRow,
                                 const ArrayView<LocalOrdinal> &Indices,
                                 const ArrayView<Scalar> &Values,
                                 size_t &NumEntries
                                 ) const =0;
 
 
    virtual void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const =0;
 
 
    virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const =0;
 
 
    virtual void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const =0;
 
    virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const = 0;
 
    virtual void leftScale (const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) = 0;
 
    virtual void rightScale (const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) = 0;
 
 
    virtual bool haveGlobalConstants() const = 0;
 
 
 
    //{@
 
    virtual const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const = 0;
 
    // TODO: first argument of doImport/doExport should be a Xpetra::DistObject
 
    virtual void doImport(const Matrix &source,
                          const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM) = 0;
 
    virtual void doExport(const Matrix &dest,
                          const Import< LocalOrdinal, GlobalOrdinal, Node >& importer, CombineMode CM) = 0;
 
    virtual void doImport(const Matrix &source,
                          const Export< LocalOrdinal, GlobalOrdinal, Node >& exporter, CombineMode CM) = 0;
 
    virtual void doExport(const Matrix &dest,
                          const Export< LocalOrdinal, GlobalOrdinal, Node >& exporter, CombineMode CM) = 0;
 
    // @}
 
 
 
 
 
    virtual std::string description() const =0;
 
    virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const =0;
 
 
    virtual void setObjectLabel( const std::string &objectLabel ) =0;
 
    // JG: Added:
 
    virtual bool hasCrsGraph() const =0;
 
    virtual RCP<const CrsGraph> getCrsGraph() const =0;
 
    // To keep the normal virtual matrix-multivector definition of apply before overloading with the region variant
    using Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node >::apply;
 
    virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X,
                      MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y,
                      Teuchos::ETransp mode,
                      Scalar alpha,
                      Scalar beta,
                      bool sumInterfaceValues,
                      const RCP<Import<LocalOrdinal, GlobalOrdinal, Node> >& regionInterfaceImporter,
                      const Teuchos::ArrayRCP<LocalOrdinal>& regionInterfaceLIDs ) const =0;
 
    // ----------------------------------------------------------------------------------
    // "TEMPORARY" VIEW MECHANISM
    void SetFixedBlockSize(LocalOrdinal blksize, GlobalOrdinal offset=0) {
 
      TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, Exceptions::RuntimeError, "Xpetra::Matrix::SetFixedBlockSize(): operator is not filled and completed."); // TODO: do we need this? we just wanna "copy" the domain and range map
 
      std::vector<size_t> stridingInfo;
      stridingInfo.push_back(Teuchos::as<size_t>(blksize));
      LocalOrdinal stridedBlockId = -1;
 
      RCP<const Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node> > stridedRangeMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(
                                                    this->getRangeMap(),
                                                    stridingInfo,
                                                    stridedBlockId,
                                                    offset
                                                    );
      RCP<const Map> stridedDomainMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(
                                              this->getDomainMap(),
                                              stridingInfo,
                                              stridedBlockId,
                                              offset
                                              );
 
      if(IsFixedBlockSizeSet()) RemoveView("stridedMaps");
      CreateView("stridedMaps", stridedRangeMap, stridedDomainMap);
    }
 
    //==========================================================================
 
    LocalOrdinal GetFixedBlockSize() const {
      if(IsFixedBlockSizeSet()) {
        Teuchos::RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::rcp_dynamic_cast<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >(getRowMap("stridedMaps"));
        Teuchos::RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::rcp_dynamic_cast<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >(getColMap("stridedMaps"));
        TEUCHOS_TEST_FOR_EXCEPTION(rangeMap  == Teuchos::null, Exceptions::BadCast, "Xpetra::Matrix::GetFixedBlockSize(): rangeMap is not of type StridedMap");
        TEUCHOS_TEST_FOR_EXCEPTION(domainMap == Teuchos::null, Exceptions::BadCast, "Xpetra::Matrix::GetFixedBlockSize(): domainMap is not of type StridedMap");
        TEUCHOS_TEST_FOR_EXCEPTION(domainMap->getFixedBlockSize() != rangeMap->getFixedBlockSize(), Exceptions::RuntimeError, "Xpetra::Matrix::GetFixedBlockSize(): block size of rangeMap and domainMap are different.");
        return Teuchos::as<LocalOrdinal>(domainMap->getFixedBlockSize()); // TODO: why LocalOrdinal?
      } else
        //TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "Xpetra::Matrix::GetFixedBlockSize(): no strided maps available."); // TODO remove this
        return 1;
    }; //TODO: why LocalOrdinal?
 
    bool IsFixedBlockSizeSet() const {
      return IsView("stridedMaps");
    };
 
    // ----------------------------------------------------------------------------------
 
    virtual void SetMaxEigenvalueEstimate(Scalar const &sigma) {
      operatorViewTable_.get(GetCurrentViewLabel())->SetMaxEigenvalueEstimate(sigma);
    }
 
    // ----------------------------------------------------------------------------------
 
    virtual Scalar GetMaxEigenvalueEstimate() const {
      return operatorViewTable_.get(GetCurrentViewLabel())->GetMaxEigenvalueEstimate();
    }
 
    // ----------------------------------------------------------------------------------
#ifdef HAVE_XPETRA_KOKKOS_REFACTOR
#ifdef HAVE_XPETRA_TPETRA
#ifdef TPETRA_ENABLE_DEPRECATED_CODE
 
    virtual local_matrix_type getLocalMatrix () const = 0;
#endif
    virtual local_matrix_type getLocalMatrixDevice () const = 0;
    virtual typename local_matrix_type::HostMirror getLocalMatrixHost () const = 0;
#else
#ifdef __GNUC__
#warning "Xpetra Kokkos interface for CrsMatrix is enabled (HAVE_XPETRA_KOKKOS_REFACTOR) but Tpetra is disabled. The Kokkos interface needs Tpetra to be enabled, too."
#endif
#endif
#endif
    // ----------------------------------------------------------------------------------
 
    virtual void residual(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & X,
                          const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & B,
                          MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > & R) const = 0;
 
    protected:
      Teuchos::Hashtable<viewLabel_t, RCP<MatrixView> > operatorViewTable_; // hashtable storing the operator views (keys = view names, values = views).
 
      viewLabel_t defaultViewLabel_;  // label of the view associated with inital Matrix construction
      viewLabel_t currentViewLabel_;  // label of the current view
 
  }; //class Matrix
 
} //namespace Xpetra
 
#define XPETRA_MATRIX_SHORT
#endif //XPETRA_MATRIX_DECL_HPP