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;






    //TODO virtual=0 // TODO: Add default parameters ?
//     virtual void multiply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, Scalar alpha, Scalar beta) 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;

    // @}



    // TODO: describe of views can be done here

    //   /** \brief Return a simple one-line description of this object. */
    //   virtual std::string description() const =0;

    //   /** \brief Print the object with some verbosity level to an FancyOStream object. */
    //   virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const =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;

    // ----------------------------------------------------------------------------------
    // "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(IsView("stridedMaps") == true) RemoveView("stridedMaps");
      CreateView("stridedMaps", stridedRangeMap, stridedDomainMap);
    }

    //==========================================================================

    LocalOrdinal GetFixedBlockSize() const {
      if(IsView("stridedMaps")==true) {
        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?

    // ----------------------------------------------------------------------------------

    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
    virtual local_matrix_type getLocalMatrix () 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